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A COMPUTER PROGRAM TO DETERMINE 
RADIATING, NONADIABATIC, INVISCID FLOW OVER A BLUNT BODY 
BY THE METHOD OF INTEGRAL RELATIONS 

By Frances W. Taylor 
Langley Research Center 

SUMMARY 

A computer program (Langley program D1250) has been developed in support of 
the study of the radiating, nonadiabatic, inviscid flow properties (pressure, temperature, 
density, velocity, and enthalpy) around a blunt body in equilibrium air by use of a modi- 
fied method of integral relations. This report describes the digital computer program 
including the methods used in the digital approximation, flow charts, D1250 program code, 
instructions for the user, and a test case with input and output listings. 

INTRODUCTION 

A modified approach to the first-order approximation in the method of integral 
relations for the numerical calculation of the inviscid, adiabatic, subsonic flow field about 
a blunt-nose body traveling at hypersonic speeds as described in references 1 and 2 has 
been extended to include the radiating nonadiabatic case (ref. 3). Langley computer pro- 
gram D1250 was developed to make the calculations for the radiating nonadiabatic flow of 
air in chemical equilibrium. Results obtained in reference 3 agree with results from 
inverse and time -dependent techniques. The agreement indicates that this relatively 
simple method of solution provides an accurate description of the blunt -body flow field in 
the subsonic region. 

A description of the computer program is presented herein along with the methods 
used in the digital approximations, flow charts, program code, instructions for the user, 

j 

and a test case with input and output listings. 

SYMBOLS 

Aj coefficients in governing differential equations 

a ?7 


tangential velocity gradient 



Bi 


E i 


hV 


H 

h 


P 

Q 


4 R,77 

R 

r b 

Rfc 

Ri 


r b 

T 


u,v 


coefficients in governing differential equations 
coefficients in governing differential equations 
coefficients in governing differential equations 

functions appearing in governing partial differential equations which are 
differentiated with respect to y 

total enthalpy, h + — * 

static enthalpy 

functions appearing in governing partial differential equations which are 
differentiated with respect to x 

pressure 

Rn 

curvature of body, 

K b 

total radiation heat flux 

universal gas constant 

body radius of curvature at x = 0 

local body radius of curvature (see fig. 1) 

weighted heat flux 

radius measured from axis of symmetry of body (see fig. 1) 

normal distance from axis of symmetry to a point on body (see fig. 1) 

temperature 

free -stream velocity 

velocity components in x- and y-direction, respectively 
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resultant velocity 


W REF molecular weight of cold air, 28.96 g/mole 
x coordinate along body surface (see fig. 1) 

y coordinate normal to body surface (see fig. 1) 

Z compressibility factor 

z axial coordinate (see fig. 1) 

6 shock displacement distance (see fig. 1) 

6 0 value of shock displacement distance at axis of symmetry 

r / transformed y-coordinate, y/6 

6-fo body inclination angle (see fig. 1) 

p density 

c o shock -wave angle (see fig. 1) 

Subscripts: 

j used (as first subscript) to denote a particular governing differential equation: 

1 shock-geometry equation 

2 continuity equation 

3 x-momentum equation 

4 y -momentum equation 

5 energy equation 

s shock -oriented properties 

stag stagnation condition 
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V 


denotes location within shock layer: 


0 body surface 

1 at shock wave 

oo free -stream conditions 

When double subscripts are used, the first (j) denotes the governing equation num- 
ber and the second (rj) denotes location within shock layer. 

Primes denote dimensional quantities. 



PROBLEM DESCRIPTION 

The equations which govern inviscid, radiating, nonadiabatic steady flow of equi- 
librium air over a blunt body traveling at hypersonic speeds are a system of nonlinear 
partial differential equations derived from the laws of conservation of mass, momentum, 
and energy. The modified method of integral relations, as described in reference 3, is 
used to transform the governing equations into a set of ordinary differential equations 
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that are numerically integrated to yield the details of the thermodynamic and flow prop- 
erties within the shock layer. The equations appear in the form used in reference 1. A 
slightly different notation is used in reference 3. The ordinary differential equations 
which form the basis for Langley program D1250 are as follows: 


H = -Ej = (1 + Q6)tan(w - 9^ 


dw 

dx 


- A 4 S + E 4 

Ba 


«2,0 - A 2i + B 2ff + E : 

dx C 

dl 


3,0 dI 2,0 


dx 


^5,0 ~ A 5 fe + B 5 af ♦ E 5 

dx C 


(1) 


The equations for the fluxes of mass, momentum (in the x-direction), and energy at the 
body surface are 


*2,0 “ P0 U 0 

* 3 ,o = p o + p o u o 2 / ^ 


*5,0 " p o u o H o 

J 


The quantities have been nondimensionalized as follows: 


x = 


y = 


u = 


V = 


x' 

g _ _5_ 

p_ P 

r b 

r b 



p-£. 

h = — SL 

r b 

P' 

oo 

Kf 

u' 

-pfl 

II 

& 

h- h' 


r b 

Kf 

v’ 

r' 


U’oo 

r = R B 



where primes denote dimensional quantities. 
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The coefficients of the governing differential equations are developed in reference 1. 
The expressions for the coefficients are 


/ v 26 cos 0 h / 

A2 = 3(l 2;0 -I 2;1 ) + ^—fe(l 2>0 -I 2!l 


A 4 = - 3 + 


26 cos 0, 


, x 26 cos x 

a 5 - 3(15,0-15,1)+— ^(15,0-15,1) 


B< = 6 3 + 


26 cos 


r b / ^ 


(j = 2,4,5) 


5 COS 0)3 


C = 6 3 + 


E -0+ Q5) 
1 Q 


( w - 0b)! 


F _ 36 /t . T \ dr b , J o , 25 cos fibV 12 * 1 d0 b 

e 2 " r^2,0 + r 2,lJ dF + 6 l 3 + r£ dF 


6^sin 0b 


( X 2,0 + 2I 2,l)dF + + + — r^~) G 2, 1 


E 4 = ^ l4 ’ 1 ^ + T + ' 


26 cos 0b\ aI 4,l d0b 26 2 sin 0b d0b 


/ 9Sb 


dx Tv 


5 cos 0b' 


+ 6 d + 5Q)(1 + ?b y G 4, 1 ■ 6G 4 ,0 - 3(5Q + + *3,1 


6 cos 


25^Q cos 


/ x 36 cos 0u/ o 

(13,0 + 2l3,l) + — K u 0 + p l u l 


+ L^!A/p oU() 2 + 2pjUj2 


E 5 " ^5,0 + r 5,l)dF 7~h,0 + 2I 5, l) dF 

+ 5(3 + ^ + efi + ^ 2 fA)(i + «Q)G 5 ,1 + 6R1 


Ic A + 2Ir 
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where 


G 2,l = Pl v l 
G 4,0 = P 0 
G 4,1 = Pi + Pi v i 2 
°5,1 = Pl v l H l 
h,o = P 0 U 0 
r 2, 1 = Pl u l 

1 3.0 = P 0 +P 0 U 0 2 

^l^Pl + Pl u l 2 
r 4,0 = 0 

*4,1 = Pl u l v l 

1 5.0 = P 0 U 0 H 0 
r 5, l = Pi u i H i 

R 1 - ‘I + 6Q> (‘ + )«B,1 - <*R,0 

Provision has been made in the governing equations for coupled radiating flow -field 
analysis. The radiation program, RATRAP (ref. 4), is an independent program which 
was incorporated into D1250 as subprograms. 

The governing differential equations are solved by a fourth -order Runge-Kutta inte- 
gration technique (see appendix A) to give shock-layer thickness, shock angle, and the 
fluxes of mass, momentum, and energy at the body surface. The numerical integration 
proceeds from the stagnation point (x = 0); however, direct substitution of the initial val- 
ues into the governing differential equations at x = 0 results in indeterminate (0/0) 
expressions. (See ref. 3.) Application of L'Hospital's rule yields the proper starting 
values for the nonzero derivatives: 
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djA (3t 35 0 t 6 0 2 )(p stas -Pj+Tj) 
dx /x-0 «o( 3 + 26 o)( 1 - Pl) v l 

'«2,0\ 6<>(3 + M* ' P l)(*) x=0 ' ( 3 + 35 “ + ^ 

r^Lo -5 0 (3 + 8 0 ) 




'dl 2> 0\ 

' <* 4=0 


where 


H 


stag - (Hl) x=0 - 


3^1, stag 


8 0 ( 3 + 80A 


Ml, 


2>0 | 


dx 


'x=0 


_ p stag a 0, stag 

5 stag 


R 


1, stag 


t 1 + 5 o) q R,l " q R,o| 


J x=0 


and p is the input value PEEP (ref. 3, eq. (36)). These equations contain 5 0 , the shock 
displacement distance at the axis of symmetry (x = 0). Since 6 0 is initially assumed, a 
correct value has to be established for the flow -field solution. The regularity conditions 

/ dI 2,0 \ 

at the sonic point, that is, a sonic velocity Uq and a zero mass flux derivative I — = 01 
provide the convergence criteria for the proper 6 0 . ' 

K Newton-Raphson iteration technique (ref. 5) is used to obtain the set of properties 
on the body which are consistent with the integrated fluxes. Rankine-Hugoniot relations 
and the shock angle are used to compute the properties immediately behind the shock. 

The fluxes pu, puv, and puH and the resultant velocity are assumed to vary linearly 
across the shock layer. Thus, the post-shock and the body properties along with this 
assumption of linear variation within the shock layer are sufficient to generate consistent 
velocity components and enthalpies across the entire shock layer. These quantities are 
computed by the following relations: 
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V 


V 


hv 

\v 


H - l5 ’ V 

%v 





The pressure distribution is then calculated from the quadratic expression 


= p 0 + (sQpoV) 7 ? + ( p i ' p o ‘ 6 Qp o u o 2 ) 7?2 

The remaining state properties are computed from the equilibrium -air correlations of 
reference 1. 

The integration continues until the pressure at the body diverges, either increasing 
or decreasing rapidly. This divergence indicates that the sonic point on the body, which 
is a singular point, has been reached. An increasing divergence of the pressure indicates 
that the initial 6 0 was too large, and a decreasing divergence of the pressure indicates 
that 5 0 was too small. The program is stopped after a specified number of iterations 
on 6 0 or when a maximum value of x, which is an input, has been reached. 


PROGRAM ORGANIZATION 


D1250 Labeled COMMON 

The following list contains all the FORTRAN variables appearing in labeled 
COMMON, except those used by the radiation package, in the order and with the dimen- 
sional information used in program D1250: 


COMMON 

label 

FORTRAN 

variable 

Description 

/REAI/ 

IREAD 

Set to 0 in EQURATS so that the NAMELIST input in 
MAIN will only be read the first time through MAIN 

/MAINP/ 

AA1(6) 

Coefficients used in the governing differential equa- 
tions, computed in COEF 


AC 

Program variable computed in MAIN 


AO 

a 0 


9 



COMMON 

label 

FORTRAN 

variable 

Description 

/MAINP/ 

AO VERB 

Ratio of semiminor axis to semimajor axis of an 
ellipse 


A1 

a l 


BB1(6) 

Coefficients used in the governing differential equa- 
tions, computed in COEF 


CAPHO 

H 0 


CAPH1 

H 1 


CAPQ 

Q 


CC 

Program variable computed in MAIN 


CC1 

Coefficient used in the governing differential equa- 
tions, computed in COEF 


CCI 

Initial computing interval used by INTI 


CCUVAR(8) 

Double -precision array used by INTI to store values 
of independent and dependent variables at which 
derivatives will be evaluated in DERSUB 


CII 

Initial interval away from singular point of stagnation 
streamline 


CIMAX 

Absolute value of maximum computing interval used 
by INTI 


COSTB 

COS 


COSW 

COS 00 


CWMTB 

COSTCO - e^j 
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COMMON 

label 

FORTRAN 

variable 

Description 

/MAINP/ 

DAF 

6 cos Qu. 

— — — - computed in DERSUB, used in COEF 


DDELDX 

d6 

dx 


DDER(8) 

Array in which derivatives are stored 


DD1 

Coefficient used in governing differential equations, 



computed in COEF 


DELL 

Lower limit for 6 


DELS 

6 2 


DELTA 

5 


DELTN 

So 


DELU 

Upper limit for 5 


DIODX(6) 

DI1DW(6) 

dI 3,0 

dx 

ot. i (J = 2,3,4, 5) 

9I 1,1 
90) 


DOSRB 

computed in DERSUB, used in COEF 


DP1DW 

dp x 

do> 


DRBDX 

dr b 

dx 


DROIDW 

dPj 

do) 


DTBDX 

dx 


DUSDW 

dug 

do) 


DU1DW 

duj 

A, » 
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COMMON 

label 

FORTRAN 

variable 

Description 

/MAINP/ 

DVSDW 

dVg 

dco 


DV1DW 

dvj 

dco 


DWDX 

dco 

dx 


EELE1(7) 

Relative error used by INTI 


EELE2(7) 

Relative zeros, used by INTI 


EELT(3) 

Values for return to program from INTI 


EE 1(6) 

Coefficients used in governing differential equations, 



computed in COEF 


EERVAL(8) 

Array used by INTI 


EMUREF 

W REF 


EOl 

Coefficient appearing in the differential equations, 



computed in COEF 


EPSR 

Accuracy criterion for Rj iteration in MAIN 


EPS21 

Accuracy criterion for iteration on p 


GGO(6) 1 

Coefficients used in governing differential equations, 


GG1(6) J 

computed in COEF 


HONE 

H 1 


ICOSW 

Control to set COSW to exactly 0 when w = ^ 

Ci 


IERR 

Integer value supplied by INTI as an error code 
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COMMON 

label 


FORTRAN 

variable 


Description 


/MAINP/ 


IGEO 

Control in GEO for body shape: 

1 sphere 

2 ellipsoid 

3 hyperboloid 

ITEXT 

Time history print option used by INTI 

ITTEST 

Iteration limit in PROPIT 

KEKONT 

Counter for number of iterations for 6 0 

KETEST 

Number of iterations desired on 5 0 

KKKK 

Radiation option: 

0 do not compute radiation 

1 adiabatic radiation 

2 nonadiabatic radiation 
sin 6u. 

computed in DERSUB, used in COEF 

r b 

OPERA 

P 

P 

PEEP 

Coefficient in equation describing pressure variation 
in stagnation region (see ref. 3, eq. (36)) 

PIN 

Poo 

PO 

p o 

PSTAG 

Stagnation pressure 

PI 

p l 

QD 

Q6 
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COMMON 

FORTRAN 


label 

variable 

Description 

/MAINP/ 

R 

Universal gas constant 


RB 

Body radius of curvature at axis of symmetry (x = 0) 

RBX 

Body radius of curvature at any x-coordinate 

RHOIN 

Poo 

RHON 

Initial guess for p 

RHOO 

p 0 

RHOl 

p l 

R1 

Weighted heat flux 

SINTB 

sin % 

SINSQW 

o 

sin co 

SINW 

sin co 

SMALLB 

Semimajor axis of an ellipse, nondimensionalized 
with Rj} 

SMALRB 

r b 

STAENT 

Static enthalpy 

SSPEC 

Print control used by INTI (if SSPEC = 0, control 
will be returned after every acceptable integration 
step; if SSPEC * 0, control will be returned at 
specified increment of independent variable) 

SWMTB 

sin(co - 

T 

T’ 
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COMMON 

FORTRAN 


label 

variable 

Description 

/MAINP/ 

TCG 

Accuracy criterion on iteration for in RANH 


TEAN 

Initial guess for T' 


THETAB 

% 


TO 

T 0 


T1 

Tl 


UIN 

Uoo 


UINS 

u 2 

u oo 


uo 

u 0 


US 

u s 


U1 

U 1 


V 

V 


VO 

v 0 


vs 

v s 


* 

VI 

V 1 


VVAR(8) 

Double -precision array used by INTI which contains 



independent variable followed by dependent variables 


W 

O) 


X 

Coordinate along body surface 


XIO(6) 

Ij, 0 (J = 2, 3, 4, 5) 
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COMMON 

label 

FORTRAN 

variable 

Description 

/MAINP/ 

XI1(6) 



XMAX 

Maximum value of x allowed 


Z1 


Z 1 


ZO 


z 0 

/QQQ/ 

QI1 ^ 




QI2 

> 

Curve -fit coefficients to radiation distribution 


QOl 

around body, used for all but final 8 0 iteration 


Q02 

J 




QRY1 

Total radiation heat flux at shock 


QRYO 

Total radiation heat flux at body 


QSTAG 

Stagnation heat flux 


The following table is a cross-reference between labeled COMMON and D1250 sub- 
programs (where X denotes labeled COMMON appearing in each subprogram and the 
radiation labeled COMMON blocks are WON, RAD, AAA, CIONCL): 
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Subprogram 

COMMON label 

REAI 

MAINP 

QQQ 

WON 

RAD 

AAA 

CIONCL 

EQURATS 

X 







MAIN 

X 

X 

X 


X 

X 

X 

RANH 


X 






GEO 


X 






FIRST 


X 






THEP 








THER 








CONT 


X 

X 





CHER 


X 






PROPIT 

COEF 


X 






CALD 


X 






DERSUB 


X 






CHSUB 








PROD 


X 


X 

X 

X 

X 

FOFX 


X 






RAD FLUX 


X 


X 

X 

X 

X 

LINE 




X 

X 

X 

X 

CONTM 

MU 




X 

X 

X 

X 

FEMP 





X 

X 

X 

DEFIOJ 

ORDERI 

ORDERV 

INVERT 

EQUIL 





X 

X 

X 

CHGBAL 

ZHV 





X 

X 

X 

LTP 





X 

X 

X 

HAFACE 

DTLNEQ 

HTP 

TRAP1 

FLUT 





X 

X 

X 

INOCAL 

TRAP 

MASBAL 





X 

X 

X 
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D1250 Subprogram Descriptions, Flow Charts, and Listings 

Langley program D1250 is one of a series of programs developed to analyze blunt - 
body flow-field data by the method of integral relations. Since many of the same com- 
putational methods were used in the various flow -field programs, each procedure was 
written as an independent subprogram to be called as needed. The radiation subprograms 
form an independent program (ref. 4) which has been added to D1250. 

This section of the report presents a brief description of the subprograms used, a 
statement of the function of each subprogram, individual flow charts, and subprogram 
listings. Listings of the radiation subprograms are also included. A list of each sub- 
program and its description is given as follows: 

Subprogram Description 


EQURATS 

MAIN 

RANH 

GEO 

FIRST 

THEP 

THER 

CONT 

CHER 

PROPIT 

COEF 


Writes program title and calls MAIN 

Reads and writes NAMELIST input, EQRAD, initializes program 
variables, and computes stagnation -streamline solution 

Calculates post-shock conditions by use of Rankine-Hugoniot equations 

Computes body geometry 

Computes properties at first step off stagnation line 

Thermodynamics routine which solves correlation equation for 
equilibrium air 

Calculates partial derivatives of equilibrium -air properties 

Calls integration routine, performs iteration on 6 0 for sphere - 
pressure results, and writes output 

Performs an additional test on 5 0 iteration 

Iteration subroutine for body properties 

Calculates coefficients of governing differential equations 
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Subprogram 


Description 


CALD Computes governing differential equations which are used in DERSUB 

DERSUB Called by Runge-Kutta integration routine INTI to evaluate the 

derivatives 

CHSUB Called by INTI to allow the user certain logical control 

FOFX A function subprogram called by ITR1 for p Q (P s t a g) iteration 

PROD Reads and writes NAMELIST input for radiation subprograms, calls 

subprograms to compute radiation, and computes thermodynamic 
properties within shock layer 

A directed graph of the FORTRAN subprograms used in D1250 is given in the flow 

chart on the following page. 
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EQURATS 


[RADIATION' 

ROUTINES 


K 


AT 62 


TBEP 


THEP 


PROD 


MAIN 


RANH 


THER 


RANH WH FIRST 


GEO 


ITR1 M 


FOFX 


THER 


CONT 


MAIN 


PROD 


RADIATION 

ROUTINES 
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u o O O O o p-l 


EQURATS 


EQURATS writes the program identification and calls MAIN. 


PROGRAM EQURATS ( INPUT » OUTPUT *T APE5= INPUT t TAPE6=OUTPUT * PUNCH ) A 1 

A 2 
A 3 

EQUILIBRIUM WITH RADIATION PROGRAMED BY FRANCES W. TAYLOR A 4 

A 5 

COMMON /REAI/ IREAD A 6 

WRITE ( 6 t 1 > A 7 

I READ- 0 A 8 

CALL MAIN A 9 

STOP A 10 

A 11 
A 12 

FORMAT ( 1H120X58HEQUI LIBRIUM WITH RADIATION PROGRAMED BY FRANCES A 13 

1W. TAYLOR/) A 14 

END A 15- 
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MAIN 


The function of subprogram MAIN is to call the NAMELIST input, EQRAD, to ini- 
tialize the program, and to compute the stagnation -line solution. The stagnation-line 
solution is obtained by an iteration on p^ ^in MAIN Pq = P 

routine ITR1. The program is ended either after a specific number of iterations on 6 0 , 
which is controlled by the input KETEST, or when a maximum value of x has been 
reached. The flow chart for subprogram MAIN is as follows: 


stag) by usin & the librar y 
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oooonn r> o 


SUBROUTINE MAIN B 1 

DOUBLE PRECISION VVAR,CCUVAR B 2 

B 3 
B 4 

COMMON /MAINP/ AAK6) , AC, AO , AOVERB, A1 , BB1 < 6 ) , CAPHO, CAPH1, C APQ , CC , C B 5 
1C1.CCI, CCUVAR (8) ,CI I, C I MAX, COSTS, COSrf , CWMT 8, OAF , ODELD X , OOER { 8 I , DD1 B 6 
2, DELL, DELS, DELTA, DELTN, DELU , DIODX ( 6 ) , D I 1DW ( 6 ) , DOST 8, DP 1 DW , DR8DX, DR B 7 
3Q1DW,DT BOX , DU SOW, DU1DW, DVSOW, DV 1DW* DWDX.tELEl (7 I , EELE2 ( 7) , EELTI3 > , 6 8 

4EE1 (6) , ERRVAL (71 , EMUREF , EO 1 , EPSR, EPS? 1 , GGO( 6 ) , GGI ( 6 > , HONE , ICOSW, I E B 9 
5RR , I GEO, I TEXT, I TTEST,KEKON T,KE TEST, KKKK, OPERA, P, PEEP, P I N, PO , P ST AG, B ID 
6PI,QD,R,RB,RBX, RHOIN, RHON, RHOO, RHOl, Rl» SINTB,SINSQH,SINM,SMALLB,SM 8 II 
7ALRB , STAENT , S SPEC ,SWMTB ,T , TCG, T EAN, THET AB, TO , T I , J IN, J I NS, UO,US,UI, B 12 
8V , VO, VS, VI, WAR (8) ,W,X,XI0(6) ,XI1(6),XMAX,Z0,Z1 8 13 

COMMON /AAA/ A( 6, 6 ) , A12 ( 16 , 12 ) , A22 ( 7 , 7 ) , ALP ( 6 ) , ALPT( 1 6 , 5 » , BMT ( 16) , 3 14 

1C(6,16) ,CH( 16,2),CP(16) ,DEGI (6) ,H { 16 ) , J AT ( 1 6, 5 > , JPH ( 16 I , K AT ( 6 ) ,KOD 8 15 

2E( 16) ,RA( 16,2 ),RBO( 16,2), RC ( 16, 2 ) , RO( 1 6 , 2 ) , RD1 ( 16, 2 ) , RE ( 16 , 2 ) , RE 1 ( 3 16 

316,2) ,SD( 16) ,TB(3 ) ,TU( 16,2 ), TU2 (16,2) , VN( 17 ) , VNEI 16 •) , VNT( 1 7 ) , VNUC 2 B 17 
46,6 ) ,W26( 6) ,W3( 16) ,Y( 16 ),RF( 16,2) ,RC1 ( 16,2 ) ,HS, IG, IGMS, IGMSP, IGP, I B 18 
50N, IS, I SP, I SPNGP, I S PNG2 , N, NG,NP, PRESS, W27.R HO, WM,SYU, ISPNG Q 19 

COMMON /RAD/ YY ( 1 00 ) , TE E( 1 00 ) , FHV ( 2D ) , NHV, NY ,C2 , I Y B 20 

COMMON /RAD/ XNN( 7, 100) 3 21 

COMMON /RAD/ NIHVC , FHVC ( 50 ) , AHV ( 50 ) , AH VL ( 20 ) B 22 

CGMMON /RAD/ Cl ,C 3, C4,FLG, C5, FLG1 B 23 

COMMON /RAD/ YDELT , DELTD.FL 1 , FL2 B 24 

COMMON /RAD/ GEE( 8 ) , EP S ( 8 ) , NU( 20 ) , ND( 70 ) , FF ( 70) , GAMP ( 70 ) , WOL ( 2C ) , F B 25 
1HVMI20) ,FHVP(20) B 26 

COMMON /CIONCL/ F ( 100, 10) ,F2( 100, 10) ,HVL( 70) ,EP,K2,K1 , IFL, IYCON, 10 6 27 

11 ,WMI ,BIJ( 100,10) .GMINUQD, 10),GPLU( 100,10) , IAED 3 23 

COMMON /RE A I / IREAD ' 8 29 

COMMON /OQQ/ QI 1, QI2 ,Q0 1,002, QRYi , QRYO , QSTAG B 30 

DIMENSION ANS (4 ) 3 31 

EXTERNAL FOFX B 32 

NAMELIST /FQRAD/ ALT, AOVERB , CC I , C I I , C I M AX , DELTX, DUMBT I M ,F ELE 1 , EFLE 8 33 

12,EELT,EITR1,EITR2,EMUREF,EPSR,EPS21 , I ALT , l GEO, I TEX T , I TTEST ,KETEST fl 34 
2, KKKK, PEEP, PIN, QI 1, QI2, 001 ,002, R,R8, R HUM, RHOIN, R1 , SSP EC , TEAN , TCG, U '3 35 

3 1 N, XMA X B 36 

ALT ALTITUDE FOR AT62 B 37 

AOVERB A BODY GEOMETRY PARAMETER B 38 

CCI INITIAL COMPUTING INTERVAL B 39 

CIMAX ABSOLUTE VALUE OF MAXIMUM COMPUTING INTERVAL USED BY INTI B 40 

CII INITIAL LINEAR INTEGRATION STEP IN FIRST 3 41 

DELTX INITIAL GUESS USED IN ITR1 B 42 


24 



oooooooooooooooooooooonnnooo 


1 

2 


3 


EELE1 

RELATIVE ERROR 

USED 

BY 

INTI 

B 

43 

EELE2 

RELATIVE ZERO 

USED 

BY 

INTI 

3 

44 

EELT 

VALUE FOR RETURN TO PROGRAM 

USED 

BY 

INTI 

6 

45 

EITR1 

RELATIVE ERROR 

USED 

IN 

ITR1 

3 

46 

EITR2 

ABSOLUTE ERROR 

USED 

IN 

ITfU 

B 

47 

EMUREF 

MOLECULAR WEIGHT OF COLD AIR USED IN THEP 




B 

48 

EPSR 

ACCURACY FOR Rl ITERATION INI MAIN 




8 

49 

EPS21 

ACCURACY FOR RHQ ITERATION IN PROPIT 




B 

50 

I ALT 

= 1 CALL AT 62 




8 

51 

I GEO 

A TEST IN GEO FOR SPHERE , ELL I P SO I D OR HYPERBOLOID 



B 

52 

ITEXT 

TIME HISTORY PRINT OPTION =0 NO PRINTOUT 

USED 

BY 

INTI 

B 

53 

ITTEST 

ITERATION LIMIT USED IN PROPIT 




3 

54 

KETEST 

THE NUMBER OF PASSES THROUGH MAIN 




B 

55 

KKKK 

IF 0 NO RADIATION, IF 1 ADIABATIC RADIATION 

IS COMPUTED 

a 

56 


IF 2 NON ADI AB ATI C RADIATION IS COMPUTED 




B 

57 

PEEP 

PRESSURE VARIATION COEFFICIENT 




a 

58 

PIN 

FREE STREAM PRESSURE 




B 

59 

R 

GAS CONSTANT 




B 

60 

RB 

BODY RADIUS 




B 

61 

RHOIN 

FREE STREAM DENSITY 




B 

62 

RHON 

INITIAL VALUE OF RHO 




8 

63 

R1 

INITIAL GUESS FOR WEIGHTED HEAT FLUX 




8 

64 

SSPEC 

PRINT CONTROL 

USED 

BY 

INTI 

B 

65 

TCG 

ACCURACY ON ITERATION FOR DENSITY IN RANH 




a 

66 

TEAN 

INITIAL VALUE OF T 




B 

67 

U I N 

FREE STREAM VELOCITY 




8 

68 

XMAX 

WHEN X = X MA X THE PROGRAM IS STOPPED 




B 

o 

69 

XLAST=X 





D 

B 

1 o 

71 

I REA 0=1 READ+l 




8 

72 

IF (IREAD-1) 1,1,3 




a 

73 

READ (5f 

EQRAD ) 




B 

74 

IF (EOF, 

5) 10,2 




3 

75 

WRITE (6 

, EQRAD) 




B 

76 

I 11=0 





B 

77 

KEKONT = 0 





3 

78 

IF (IREAD.EQ.l) IPROD=0 




B 

79 

I RE A0= I RE AD + I 




a 

80 

R2=R1 





B 

81 

KEK0NT=KEK0NT+1 




B 

82 

IF (KEKONT. GT. KETEST) GO TO 1 




B 

83 

IF (KKKK 

. NE .2 ) RL=0 




B 

84 

IF (III. 

NE.O) GO TO 4 




8 

85 
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IF (IALT.NE.1) GO TO 4 8 86 

CALL AT 62 (ALT, ANSI B 87 

PIN=ANS(2)*478. 80258 8 88 

RH0IN=ANS(1)*. 515379 B 89 

WRITE (6,11) RHOIN, PIN B 90 

4 X=0 3 91 

W=3. 1415927/2. B 92 

THETAB=3. 1415927/2. 3 93 

3 94 

ICOSW IS 0 WHEN COS(W)=0 IC0SW=1 COS(W) IS COMPUTED B 95 

B 96 

IC0SW=0 B 97 

UINS=UIN**2 B 98 

DRBDX=1. B 99 

DTBDX=-1. B ICO 

SMALRB=0 B 101 

RBX=1. B 102 

CAPQ=1. B 103 

T=TEAN B 104 

. CALL RANH B 105 

IF (III .NE.OI GO TO 5 B 106 

R078=DUMBTIM/RH01 B 107 

DELU=1.1*R078 B 108 

DELL=.8*R078 B 109 

5 111=111+1 8 110 

DELTN=(DELL+DELU>/2. B 111 

DELTA=DELTN B 112 

QD=CAPQ*DELTA B 113 

WRITE (6,12) DELTN, XLAST 8 114 

T WO DEL=2.* DELTA B 115 

THRDEL=3.*DELTA B 116 

DELS=DELT A**2 B 117 

AC=(3.+THRDEL+DELS)/V1 B 118 

BC= DELTA* ( 3 .+TWODEL I B 119 

CC=DELTA*(3.+DELTA) B 120 

U0=0 B 121 

V0=0 B 122 

S INTB= 1 * B 123 

DIM=2.*DELTA*SINTB/DRBDX B 124 

6 P0=P1+.5*RH01*V1**2 B 125 

R1=R2 B 126 

CALL I TR1 (PO,DELTX,FOFX,EITR1,EITR2,1QO,ICODE) B 127 

IF ( ICODE.GT .0 1 WRITE (6,13) ICODE B 128 
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P1MPQ=P1-P0 B 129 

RHOO=-( AC**2*P1MP0**2 ) / (PEEP*CC*#2*P0 ) B 130 

CAPHQ=C APH1— { 3 .*R 1 ) / ( AC*P1MP0 ) B 131 

RH00A0=(AC*P1MP0)/CC B 132 

DWDX=AC*( P0— P 1+V1 ) / ( BC* ( 1 ,-RHOl ) ) B 133 

A1=C0WDX*(1 RH01 ) + 1 . ) /RHO 1 B 134 

A0=RH00A0/RH00 B 135 

DELTD=DELTA*RB B 136 

UST AR=0 B 137 

IF (KEKONT.EQ.KETEST) WRITE (6,14) UO, RHOO, AO, RH00A0 , C APHU , PO , OWD X 8 138 
1,DELTD, DELTA B 139 

I PROD= I PROD+1 B 140 

B 141 

IF KKKK=0 THE RADIATION IS NOT COMPUTED, KKKK= 1 ADIABATIC S 142 

RADIATION IS COMPUTED, KKKK=2 FOR NONADI ABAT I C RADIATION B 143 

B 144 

CALL PROD (USTAR, IPROD, AO, A1,QRY0,QRY1 ) B 145 

IF (KKKK.NE.2) GO TO 9 B 146 

R1P=(1.+DELTA)**2*QRY1-QRY0 B 147 

ABR=ABS((R1P-R1)/R1P) a 148 

IF ( ABR-EPSR) 8,8,7 B 149 

R2=R1P 8 150 

GO TO 6 B 151 

B 152 

KEKONT COUNTS THE NUMBER DF PASSES THRU MAIN B 153 

B 154 
B 155 

TEST FOR NEXT TO LAST TIME B 156 

B 157 

IF ( ( KE TEST- KEKONT) .NE. 1 ) GO TO 9 B 158 

QSTAG=1.E— 7*QRY0*RH0IN*UINS*UIN B 159 

DELTA= ( DELU+DELL ) /2 , B 160 

PST AG=PQ B 161 

CALL FIRST B 162 

0 STOP B 163 

B 164 
8 165 

1 FORMAT ( 1H06X4HAT 62 6X5HRH0 INE15. 8*6X3HPINE15.8) B 166 

12 FORMAT ( 1H06X5HDELTNE1 5.8,6X1HXE15.8) 8 167 

13 FORMAT (1 HQ6X6HIC QOE= 14,10X2 7H ITERATION DID NOT CONVERGE) B 168 

14 FORMAT ( 1 HI 5X27HST AGNAT ION POINT COND I T I0NS6X2HU0E1 5 . 8, 4X4HRH00E 1 5 8 169 

1.8, 6X2H AOE1.5. 8/35X6HRHOOAOE 15. 8 ,6X2HH0E15.8,6X2HP0E15. 8 ,4X4HDWDXE1 B 170 

25.8/6X24HSHOCK STAND-OFF D I STANCE5X5HDELTDE 15 . 8 , 1 X8HDELTA/RBE15. 8/ B 171 

3) B 172 

END B 173- 
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RANH 


RANH computes the post -shock conditions (properties, fluxes, and partial deriva- 
tives with respect to oo) with the use of the Rankine-Hugoniot equations. It is called by 
MAIN for the stagnation -line solution and then by DERSUB in each integration step. The 
flow chart for subprogram RANH is as follows: 
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SUBROUTINE RANH C I 

DOUBLE PRECISION VVAR»CCUVAR C 2 

C 3 
C 4 

COMMON /MAINP/ AAI ( 6 ) , AC, AO , AOVERB, A 1, BB1 C 6 ) , CAPHO, C APH I, C APQ , CC , C C 5 
1C1,CCI ,CCUVARC8),CI I , C I MAX , COSTB , CQSW , CWMT B, DAF , OOELDX, ODER ( 8 ) , DD1 C 6 
2, DELL, DELS, DELTA, DELTN,DELU,0I0DX{6) ,DILDW(6) ,DOSR8, DPiDW,ORBOX,QR C 7 
30 IDWtDTBDX »DUS'DW» DU1DW , DVSDW,0V1DW»DWDX,EELE1(7)» EELE2 < 7) , EELTI 3) , C 8 
4EEI (6 ), ERRVALI7), EMUREF , EO 1 , EPSR, EPS2 1 , GGO ( 6 ) , GG1 ( 6 J , HONE , ICOSW,IE C 9 
5RR, IGEO,ITEXT,ITTEST,KEKONT,KETEST,KKKK,OPERA,P,PEEP, PIN, PO, PSTAG, C 10 
6P1, QD,R,RB»RBX,RHOIN,RHGN, RHOO, RH01 , R 1 , S INTB, SINS3W,SINW, SMALLB , SM C 11 
7ALRB,STAENT,SSPEC,SWMTB,T,TCG,TEAN,THETAB,T0,T1,UIN,UINS, U0,US,U1 , C 12 
8V,V0,VS,V1,VVAR(8 >,W,X,XI0<6) , X 1 1 1 6 ) , XMAX, ZO, Z 1 C 13 

C 14 

RANKINE HUGONIOT RELATIONS FOR EQUILIBRIUM COMPOSITION C 15 

C 16 

I RHO=0 C 17 

DTBDW=0 C 18 

SINW=S INI W ) C 19 

IF (ICOSWI 1,1,2 C 20 

1 C0SW=0 C 21 

SWMTB=0 C 22 

GO TO 3 C 23 

2 COSW=COS(W) C 24 

$WMTB= S IN ( W-THETAB) C 25 

3 CWMTB=COS l W-THETAB) C 26 

SWS=SINW**2 C 27 

US=COSW C 28 

RHO 1=RH0N C 29 

4 VS=-SINW/RH01 C 30 

CAPH1=.5+3.5*PIN/(RH0IN*UINS) C 31 

HON E=C APH1 C 32 

STAENT=CAPHl-lUS**2+VS**2)/2. C 33 

CALL THEP (R,P1,RH01,T1,Z1,STAENT,EMUREF,PIN,RH0IN,UINS,KEZI C 34 

RH01P=(PIN/(RHOIN*UINS)+SWS-P1)/VS**2 C 35 

FRH01=RH01-RH01P C 36 

DELR01=.l C 37 

RH01D=RH01+DELR01 C 38 

VP=-S I NW/ RHOl D C 39 

STAP=CAPHl-(US**2+VP**2)/2. C 40 

CALL THEP (R, PD, RH01D, TX, ZX , STAP, EMUREF, PIN, RHOIN,UINS,KEZ) C 41 

RH01DP=MPIN/(RH0IN*UINS)+SWS-PD)/VP**2 C 42 


29 



co n o 


FPRH01= ( RHG1D— RH01 DP-FRH01 ) /DELR01 C 43 

AFMFP= A8S ( FRH01/FPRH01 ) C 44 

IF ( AFMFP-T CG ) 6,6,5 C 45 

5 RHO 1=RH01— FRHOI/FPRHOI C 46 

IF ( IRHO.GT . I TTEST ) GO TO 7 C 47 

IRHO= IRHO+l C 48 

GO TO 4 C 49 

6 U1=US*CWMTB-VS*SWMTB C 50 

V1=US*SWMTB+VS*CWMTB C 51 

XI 1 ( 2 ) =RH01*U l C 52 

XII (3)=P1+RH01*UI**2 C 53 

XII (4>=RH01*U1*V1 C 54 

XU (5)=RH01*U1*CAPHI C 55 

XI1(1)=0 C 56 

CALL THER IR, PI, RHO 1,T 1, 21 , STAENT , EMUREF, PIN, RHOIN, 01 NS ,CiP,C2P,Cl C 57 

1D,C2D) C 58 

DR01DW=SINW*C0SW*(C1D*< l.-l./RH01**2)-2.*(l.-l./RH01) ) / ( ( S WS/RHOl * C 59 

1*2 ) *( 1 •— C 1D/RH01 ) — C20 ) C 60 

DSH1DW=SINW*C0SW*(1.-1./RH01**2)+SWS/RH01**3*DR01DW C 61 

0P1DW=C1D*DSH10W+C2D*DR01DW C 62 

DUSDW=-S I NW C 63 

DVSDW=- 1 . /RH01*I US+VS/RH01 *DR010W ) C 64 

DU1DW=0US0W*CWMTB— DVSOW*SWMTB-US*SWMTB*{l.-OTBDW)-yS*CWMTB*( 1.-DT8 C 65 

1DW) C 66 

DV1 DW=OUS DW*S WMTB+OVSOW*CWMT B-HJS*CWMTB*(1 . — DTBDW)— VS*SWMTB*( 1 .-DTB C 67 

1DW) C 68 

DI 1DW( 2 ) =RH01 *0U1 0W+U1*0R0 1 DW C 69 

DI 1DW( 3 ) =DP 1DW+2. *RH01*U1*0U10W+U1**2*0R010W C 70 

DI10WI4)=RH01*V1*OU1DW+RH01*U1*DV10W+U1*V1*OR01DH C 71 

DI1DW( 5)=CAPH1*0I 10W(2 » C 72 

RETURN C 73 

7 WRITE (6,8) C 74 

STOP C 75 

C 76 
C 77 

FORMAT ( 1 HG6X39HITERAT ION LIMIT IN RANH PROGRAM STOPPED) C 78 

END C 79- 
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GEO 


Subprogram GEO computes the body geometry. The user has a choice of three 
body shapes which are controlled by an input IGEO. IGEO = 1 computes sphere geometry, 
IGEO = 2 computes ellipsoid geometry, and IGEO = 3 computes hyperboloid geometry. 
GEO is called by FIRST for the stagnation-line solution and by DERSUB for each integra- 
tion step. The flow chart for subprogram GEO is as follows: 
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SUBROUTINE GEO F 1 

C F 2 

DOUBLE PRECISION VVAR,CCUVAR F 3 

C F 4 

COMMON /MAINP/ AA 1 ( 6) , AC, AO , AOVERB , A1 , BB1 ( 6 ) , CAPHO, CAPH1 , CAPQ,CC , C F 5 
ICI, CCI ,CCUVAR(8),CII,CIMAX,C0STB,CQSrt,CWMTB,DAF,DDELDX,DDER(8),DDl F 6 
2 » DELL , DELS * DELT A, DELTN , DELU , DIODX I 6) , D I 1DW ( 6) , OOSRB , DP IDW, DRBDX, DR F 7 
301DW,DTBDX,DUSDW,DUlDW,DVSDW,DVlDW,D?iDX,EELEim,EELE2(7),EELT(3) , F 8 
4EE1 (6) , ERRVAL (7), EMUREF ,E0 1 , EPSR, EPS2 1 , GGOI 6 ) ,GG1 ( 6 ) , HONE , ICOSW,IE F 9 
5RR, IGEO»ITEXT,ITTEST,KEKONT,KETEST,K<KK., OPERA, P, PEEP, PIN»PO,PSTAG, F 10 
6P1,QD,R,RB,RBX,RH0IN,RH0N, RHOO, RH01,R1 , SINTB, S I NSQW , S INirf, SMALLB, SM F II 
7ALRB,STAENT»SSPEC,SWMTB,T,TCG,TEAN,THETAB,T0,T1,UIN,UINS,U0,US,U1, F 12 
8V,V0,VS,V1,VVAR(8) , W» X , XI 0 { 6 > , X II { 6 ) , XMAX , ZO, Z1 F 13 

F 14 

IGEO=l COMPUTE SPHERE GEOMETRY IGE0=2 COMPUTE ELLIPSOID GEOMETR F 15 
I GE0=3 COMPUTE HYPERBOLOID GEOMETRY F 16 

F 17 

SMALLB= AOVERB F 18 

IF (IGEO-2I 1,2,3 F 19 

F 20 

SPHERE - CONE F 21 

F 22 
F 23 

SMALRB= S I N C X ) F 24 

THETAB=3. 14159265/2. -X F 25 

DRBDX=COS ( X ) F 26 

OTBOX=-l. F 27 

RBX=1. F 28 

GO TO 4 F 29 

ELL IPSIOD F 30 

BSQ=SMALLB**2 F 31 

RSQ=SMALRB**2 F 32 

BSMRS=BSQ-RS9 F 33 

RACE=RSQ/ BSMRS F 34 

ZB=AOVERB*SMALLB-AOVERB*SQRT(BSMRS) F 35 

AAA=AOVERB*SMALLB F 36 

PTAN={(SMALLB/AAA)* I AAA-ZB ))/ISQRT(2. *AAA*ZB-ZB**2 ) ) F 37 

THE TAB = AT AM PT AN) F 38 

RBX= ( ( 1 .+AOVERB**2*RAC E ) *# ( 3 . /2. ) ) / ( AOVERB/SQRTI BSMRS I * ( 1 . + RACE ) ) F 39 
CAPQ=1./RBX F 40 

DTBDX=-CAPQ F 41 

GO TO 4 F 42 
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HYPERBOLOID F 43 

F 44 

BSQ=SMALLB**2 F 45 

RSQ=SMALRB**2 F 46 

BSPRS=BSO+RSQ F 47 

RACH=R SQ/ 8SPRS F 48 

AAA=AOVERB*SMALLB F 49 

SQBSPRS=SQRT ( BSPRS) F 50 

ZB=-AAA+AOVERB*SQBSPRS F 51 

PTAN=SQB$PRS/ (A0VERB*SMALRB) F 52 

THETAB=ATAN(PTAN) F 53 

RBX=((1.*A0VERB**2*RACH)**( 3./2. ))/( AOVERB/SQBSPRS* { 1 .-RACHH F 54 

4 CAPQ= 1 . /RBX F 55 

QO=CAPQ*DELT A F 56 

RETURN F 57 

END F 58- 
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FIRST 


Subprogram FIRST is used to compute the properties at the first step off the stag- 
nation line. The integration cannot start at the stagnation point since it is a singular 
point. The flow chart for subprogram FIRST is as follows: 



Compute properties 
for first step off 
stagnation streamline 

x > Uq* H 0> r b 
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SUBROUTINE FIRST 1- * 

DOUBLE PRECISION VVAR,CCUVAR L 2 

L 3 
L 4 

COMMON /MAINP/ AA1 ( 6) , AC, AO , AOVERB, A1 , BB1 ( 6 ) , CAPHQ, CAPH1, C APQ, CC , C L 5 

1C1,CCI , CCUVAR(8),CI I ,C I MAX, COSTS, COStf , CWMT B, DAF, DDELDX , DOER 18> ,DD1 L 6 

2, DELL , DELS » DELTA, DELTN, DELU , DI0DXI6),D11DW(6) ,DOSRB , DP1DW, DRBOX, DR L 7 

301DW,DTB0X,DUSDW,DUlDW,DVSDW,DVlDW,DtfDX,EELEim,EELE2m, EELTI3I , L 8 

4EE1 (6) , ERRVALC7), EMUREF, EOI, EPSR, EPS2 1 , GGO ( 6) , GG1 ( 6 ) , HONE, ICOSW.IE L 9 

5RR»IGEO»ITEXT,ITTEST tKEKONT * KETEST ,KKKK,OPERA,P ,PEEP, PIN, PO,PSTAG , L 10 

6P1, QD, R,RB,R8X,RHQIN,RH0N, RHOO, RH01, Rl, SINTB, SINSQW , S INW, SMALLB , SM L 11 

7ALR B» STAENT,SSPEC,SWMTB,T,TCG,TEAN,THETAB,T0,T1,JIN,UINS, UO ,US ,U1 , L 12 

8V»V0,VS,V1»VVAR(8>,W,X,XI0(6),XI1(6),XMAX,Z0,Z1 L 13 

C L 14 

X=CII L 15 

IC0SW=1 L 16 

U0= AO*C II L 17 

W=W+DWDX*C 1 1 L 18 

CAPHO=CAPHO+UO**2/2. L 19 

SMALRB=CI I L 20 

CALL GEO L 21 

THETAB=3. 14159265/2. +DTBDX*CII L 22 

DELTA=DELT A— ( Cl .+QD )/C APQ ) * ( W-T HETAB ) *DTBDX*C 1 1 L 23 

XIO( 3 ) =P0+RH00*U0**2 *- 24 

X 1 0 C 4 1 = RHOO*VO*UO <- 25 

XIQ(5)=RH00*U0*CAPH0 L 26 

XIO<2)=RHOO*UO L 27 

1 CALL RANH L 28 

CALL CONT L 29 

RETURN >- 30 

END L 31- 
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THEP 


THEP is the thermodynamics routine that solves the correlation equations for equi- 
librium air (ref. 1). THEP is called by RANH, PROPIT, and FOFX. It has a calling 
sequence in which R, EMUREF, and the free -stream conditions PIN, RHOIN, and UINS 
are constants. The values of density RHODUM and static enthalpy STAENT are given at 
the particular point at which THEP is called. The corresponding pressure PDUM, tem- 
perature TDUM, and compressibility ZDUM are calculated. There is an indicator KEZ 
which is set to 1 if ”RHO IS SMALLER THAN ALLOWED.” This is only used in PROPIT. 
When KEZ = 1, PRO PIT returns immediately to CHER. CHER then sets DELL = DELTN 
and calls MAIN. The program starts again with a new 5 0 . The flow chart for sub- 
program THEP is as follows: 
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SUBROUTINE THEP ( R* PDUM* RHODUM , TDUM f ZDUM* 

STAENT tEMUREFf PINf RHOIN^U 

K 

1 

1 1 NS f KEZ ) 


K 

2 



K 

3 



K 

4 

KEZ WAS ADDED TO PROGRAM 12-14-67 


K 

5 

IF RHO IS SMALLER THAN ALLOWED OR HH IS 

0 OR NEG* WHEN THEP IS 

K 

6 

CALLED FROM PROPIT KEZ^I AND THE PROGRAM 

RETURNS TO CHER 

K. 

7 

IN CHER DELL=DELTN MAIN IS CALLED 


K 

8 



K 

9 



K 

10 

KEZ— 0 


K 

11 

HUM=EMUREF*UI NS*STAENT/R 


K 

12 

RHO^RHOIN/ 1* 29313E— 3*RHGDUM 


K 

13 

IF <RH0*LE*.00G03) GO TO 11 


K 

14 

IF ( RHO *GE *30 • ) GO TO 13 


K 

15 

IF <HUM*LE.O. I GO TO 10 


K 

16 

IF { HUM*GT .5800* ) GO TO 1 


K 

17 

P=( *975134E— 3 >*HUM*RHO 


K 

18 

Z = 1*0 


K 

19 

GO TO 9 


K 

20 

IF (HUM*GT *10500® ) GO TO 2 


K 

21 

P=*00345*RHG*EXP( *8 54* ALOG C HUM) I 


K 

22 

Z=1.0 


K 

23 

GO TO 9 


K 

24 

IF ( HUM *GE*35500* ) GO TO 3 


K 

25 

RHOLOG= AL0G10 (RHO ) 


K 

26 

XD=5**AL0G10( HUM! —20* 


K 

27 

PLOG=* 955+RHOLOG* I * 1545 + *0 1 3 1*RH0LQG ) *XD+ 

*016*RHQL0G*XQ*( 2*75-XD) 

K 

28 

P=EXP( 2»302585*PLOG J 


K 

29 

GO TO 5 


K 

30 

DV= 1*56 5+1* 036* ALOG 10 (RHO ) 


K 

31 

IF (HUM *GE • 178000 * ) GO TO 4 


K 

32 

XD=5.*ALGG10(HUM)-24* 


K 

33 

PLGG=*1336*XD+*00934*XD**3+DV 


K 

34 

P^EXP ( 2 *302 5 85* PL OG ) 


K 

35 

GO TO 5 


K 

36 

IF ( HUM *GE*600000 • ) WRITE (6,17) 


K 

37 

PLOG= * 95* ALOG 10 ( HUM 1-4 • 58+DV 


K 

38 

P=£XP(2*302585*PLOG) 


K 

39 

GO TO 5 


K 

40 

COMPUTE Z(H,RH0) 


K 

41 

X 1= 1 *0- ALOG 10 ( RHO ) 


K 

42 
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XD= ALQG10 « HUM )-4. +*052*XI- .004*XI **2 

K 

43 


IF ( XO • GT «0 • ) GO TO 6 

K 

44 


Z = 1.0 

K 

45 


GO TO 9 

K 

46 

6 

IF (XD,GT.0.55) GO TO 7 

K 

47 


Z=1 .0+.53*XD**2 

K 

48 


GO TO 9 

K 

49 

7 

IF (XD.GE.1.3) GO TO 8 

K 

50 


X 1= 1* 3— XD 

K 

51 


Z = 2<.0 + XI*<-1*78+XI*(.21+XI*( 1.09-*446*XI*XI 1 ) ) 

K 

52 


GO TO 9 

K 

53 

8 

IF (XD.GE.2.) WRITE 16,17) 

K 

54 


X 1 = 1 * 9-XQ 

K 

55 


Z = 3.831+XI*(-5.G19+XI*( 3.41+.24*XI ) ) 

K 

56 

9 

T=(273.*P)/(Z*RHG I 

K 

57 


ZOUM=Z 

K 

58 


TDUM=T 

K 

59 


PDUM=1.G1325E6*P/(RH0IN*UINS) 

K 

60 


GO TO 14 

K 

61 

10 

WRITE (6,16) 

K 

62 


GO TO 12 

K 

63 

11 

WRITE (6,18) 

K 

64 

12 

KEZ=1 

K 

65 


GO TO 14 

K 

66 

13 

WRITE (6,19) 

K 

67 


GO TO 15 

K 

68 

14 

RETURN 

K 

69 

15 

STOP 

K 

70 

C 


K 

71 

C 


K 

72 

16 

FORMAT U0X22HH IS ZERO OR NEGATIVE ) 

K 

73 

17 

FORMAT ( 10X25 HH IS LARGER THAN ALLOWED ) 

K 

74 

18 

FORMAT ( 10X28HRH0 IS SMALLER THAN ALLOWED ) 

K 

75 

19 

FORMAT (10X27HRH0 IS LARGER THAN ALLOWED ) 

K 

76 


END 

K 

77- 
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THER 


Subprogram THER is called by RANH to calculate the partial derivatives of 
equilibrium -air properties which are used in computing the derivatives at the shock 
wave. THER has a calling sequence in which R, EMUREF, and the free-stream condi- 
tions PIN, RHOIN, and UINS are constants. PDUM (PI), RHODUM (RHOl), and STAENT 
are computed in RANH. The quantities C1P (dp/dh) and C2P (dp/d p) are calculated in 
THER. The flow chart for subprogram THER is as follows: 



39 







o o o o 


SUBROUTINE THER C R, POUM *RHOOUM» TDUM » ZDUM, STAENT f EMUREF, PIN,RHOIN,U I 1 

HNSfClPfC2P»CiD*C201 I 2 

I 3 

THIS IS THE NEW THERMOCHEMICAL ROUTINE THER 4-19-67 I 4 

I 5 
I 6 

Z = ZDUM I 7 

P~P OUM I 8 

HUM=EMUREF*UI NS*STAENT/R I 9 

RHO=RHGIN/l*29313E-3*RHODUM I 10 

RUN=1. I 11 

XC1 PDP-5* *ALOG10( HUM) -20* I 12 

IF ( RHO.L E* *00003 ) GO TO 8 I 13 

IF (RHO *GE. 30* ) GO TO 9 I 14 

IF (HUM*LE.G.J GO TO 6 I 15 

IF ( HUM *GT • 5800* ) GO TO 1 I 16 

ClP=.975134E-3/RUN*RH0 I 17 

C2P=*975134E-3*HUM I 18 

GO TO 5 I 19 

1 IF (HUM. GT. 10500. I GO TO 2 I 20 

C1P=.345E“2*.854*RH0/(RUN*HUM**.146) I 21 

C2P=.345E-2*HUM I 22 

GO TO 5 I 23 

2 IF (HUM *GE *35500* ) GO TO 3 I 24 

ClP=5.*P/(HUM*RUN)*(.1545+( *0131+. 016^(2*75-2. *XC1PDP n*ALOG10(RHO I 25 

in i 26 

C2P=P/RH0*(1.+XC1PDP*( *0131+.016*(2.75-XClPDP)n I 27 

GO TO 5 I 28 

3 DV=1.565+1.G36*ALGG10(RH0) I 29 

IF (HUM.GE. 178000.) GO TO 4 I 30 

C2P=1.036*P/RHO I 31 

ClP=5.*P/{HUM*RUN)*(.1336+.02802*(5.*ALQG10{HUM)-24. >**2) I 32 

GO TO 5 I 33 

4 IF (HUM. GE. 600000* ) GO TO 7 I 34 

C1P=.95*P/(HUM*RUN) I 35 

C2P=1.036*P/RH0 I 36 

5 T=( 273 * *P ) / ( Z*RHO ) I 37 

T OUM=T I 38 

C10=C1P*1 .01325E6/ ( RHO I N*R ) *EMUR£F I 39 

C2D=C2P*1 *01 325E6/I 1. 29313E-3*UI NS ) I 40 

GO TO 11 I 41 

6 WRITE ( 6 t 12 ) I 42 
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GO TO 10 I 43 

7 WRITE <6,13) I 44 

GO TO 10 I 45 

8 WRITE (6,14) I 46 

GO TO 10 I 47 

9 WRITE (6,15) I 48 

GO TO 10 I 49 

10 CONTINUE I 50 

11 RETURN I 51 

I 52 

I 53 

12 FORMAT (10X22HH IS ZERO OR NEGATIVE ) I 54 

13 FORMAT (10X25HH IS LARGER THAN ALLOWED ) I 55 

14 FORMAT ( 10X28HRH0 IS SMALLER THAN ALLOWED ) I 56 

15 FORMAT (10X27HRH0 IS LARGER THAN ALLOWED ) I 57 

END I 58- 
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CONT 


In subprogram CONT the array of variables for INTI, the integration routine, is 
established and INTI is called to evaluate the derivatives. CONT writes the output from 
the integration routine. PROD is called to compute the 77 distribution and if KKKK is 
greater than zero, to compute the radiation. To minimize the use of the time-consuming 
radiation subroutine RAD FLUX, curve fits to the radiation distribution around the body 
are used for all except the last iteration. The curve-fit expressions are 

QRY1 = (1 + QI 1 x X 2 + QI 2 X X 4 )QRYl stag 

QRYO = (1 + QOl X X 2 + Q02 x X^QRYOgtag 

where Qll, QI2, QOl, and Q02 are coefficients which describe the curve fit. The curve- 
fit coefficients are input values and are based on available data for radiation distributions 
or on the experience of the user. If no basis exists for these coefficients, the radiation 
heat fluxes should be computed at each integration step. The flow chart for subprogram 
CONT is as follows: 
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SUBROUTINE CONT G I 

DOUBLE PRECISION VVAR,CCUVAR G 2 

G 3 
G 4 

COMMON /MAI NP/ AA 1 { 6 I , AC? AO , AOVERB » A1 f BB1 ( 6 ) , CAPHO , C APH1 * C APQ, CC , C G 5 
1C 1 , CCI f CCU VAR *8),CIX,CIMAX,COSTB f COSW,CWMTB t DAF f DDELD X* DOER < 8 ) ? DO 1 G 6 
2,DELL,DELS,OELTA,DELTN?DELU,DIOQX{6) 9 DIIDWC 6) , DOSRB, DP 1 DW , DRBDX, OR G 7 
301DW, DT8DX, DUSDW, OU1DW , DVSQW»DV1DW,DWDX,EELE1(7),EELE2{7) , EELTC 3 ) , G 8 
4EE1 (6) ,ERRVAL(7) f EMUREF , EOl , EPSR f EPS2 1 , GGO ( 6 > , GG1 { 6 ) , HONE * ICGSW, IE G 9 
5RR, I GEO, I TEXT, I TT EST , KEKQN T , KETES T , < < KK , OPERA , P,PEEP, PI N, PO , PSTAG , G 10 
6P1, QOt RtRBtRBXfRHOINf RHONt RHOOt RHOltRl , SINTB, S I NSQW , S INW, S MAL LB , SM G 11 
7ALRB v STAENT»SSPECf SWMTByTyTCGvTEANyTHETABf T0*T1 ,UIN , U I NS, UO ,US ,U1 , G 12 
8V,V0,VS,V1,VVARC8I ,W,X, XI0C61,XI1 (6), XMAX, Z0,Z1 G 13 

COMMON /QQQ/ 01 1 , QI 2 9 QO 1 , Q02 , QRY1 , QRYO , QSTAG G 14 

EXTERNAL DERSUB G 15 

EXTERNAL CHSU8 G 16 

PCOREL= 1 « G 17 

I T=0 G 18 

WAR C 1 I = X G 19 

WAR ( 2 ) =DELTA G 20 

WAR ( 3 1 =W G 21 

DI ODX ( 4 1=0 G 22 

DO 1 KM=2,5 G 23 

VVAR(KM+2)=XI0CKM) G 24 

CONTINUE G 25 

IF I GEO IS EQUAL TO 2 DRBOX IS INTEGRATED G 26 

NE IS THE NUMBER OF EQUATIONS TO BE INTEGRATED G 27 

G 28 

IF HGEG-2) 2,3,3 G 29 

2 N E= 6 G 30 

GO TO 4 G 31 

3 N E= 7 G 32 

VVAR(85=SMALRB G 33 

4 IF ( KEKONT *LT*KETEST1 GO TO 5 G 34 

WRITE (6,11) (VVARC JI,J=1,8) G 35 

5 CALL INTI < IT , NE, 3, CCX , SSP EC , C I MAX , I E RR » WAR , CCUV AR , DDER, E E LEI , EE L G 36 

1E2,EELT,ERRVAL,DERSUB,CH$UB,ITEXT) G 37 

COSTB=COS (THETAB) G 38 

SINTB=SIN(THETAB) G 39 

SMALLRS = ( 1 ® -e-DELTA I^COSTB G 40 

ZS=1*-(1.*DELTA)*SINTB G 41 

DELTD=DELTA*RB G 42 
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I PR 00= 2 G 43 

UST AR=2 • G 44 

C G 45 

C IF KKKK=0 THE RADIATION IS NOT COMPUTED, KKKK=1 ADIABATIC G 46 

C RADIATION IS COMPUTED, KKKK=2 FOR NONADI ABATIC RADIATION G 47 

C G 48 

IF (KEKONT.LT. KETEST) GO TO 6 G 49 

WRITE (6,12) ( CCU VAR ( JK)»JK=1,8), (XU (JJ)»JJ=2,5) » R8X , THET AB ,C APQ , G 50 
1(D0ER( ID ,IL=2,8) G 51 

CALL PROD (USTAR, IPROD, U0,U1 , QRYO, QRY l ) G 52 

IF (KKKK.NE.2) GO TO 7 G 53 

R1=(1.+QD)*(1.+DAF) *QRY 1-QRYO G 54 

QRY0P=1 .E-7*QRY0*RH01N*UINS*UIN G 55 

QRATIO=QRYOP/QSTAG G 56 

GO TO 7 G 57 

> XS=X**2 G 58 

QST AR1 = ( 1 • +Q 1 1*XS +Q I2*XS**2 ) *QRY1 G 59 

QSTAR0=(l.*Q01*XS+Q02*XS**2)*QRY0 G 60 

IF (KKKK.NE.2) GO TO 7 G 61 

Ri=(l.+QD)*(1.+DAF)*QSTAR1-QSTAR0 G 62 

7 PRATIO=PO/PST AG G 63 

PDI F=PCOREL— PRATI 0 G 64 

IF (KEKONT.EQ. KETEST) WRITE (6,13) PRATIO, QRATI 0 G 65 

IF (PDIF) 8,9,9 G 66 

8 DELU=DELTN G 67 

CALL MAIN G 68 

9 PCOREL=PRATIO G 69 

IF (X-XMAX) 5,5,10 G 70 

10 RETURN G 71 

C G 72 

C G 73 

11 FORMAT ( 1H01 1X44HPR0PERTIES AT FIRST STEP OFF STAGNATION L INE/8X1H G 74 

1XD1 5. 8 , 5X5HDELT ADI 5 .8,9X1HWD15. 8, 5X5HIO(2)D15.8, 5X5H 10 ( 3)015.8/5X5 G 75 
2HI0(4)D15.8,5X5HI0(5)D15.8, 4X6HSMALRB0 15. 8/ ) G 76 

12 FORMAT (1H1/10X25HBODY AND SHOCK PR0PERTIES/9X1HXD15. 8, 5X5HDELTAD1 G 77 

15.8,9X1HWD15.8,5X5HI0( 2)01 5.8, 5X5HI 0(3) D15.8/5X5HI0 (4) 015.8, 5X5HI0 G 78 
2(5)015. 8, 4X6HSMALRBD15 . 8, 5X5HI1(2)E15.8,5X5HI1(3)E15.8/5X5HI1(4)E1 G 79 
35.8»5X5HI1(5)E15.8» 7X3HRBX E15 .8»4X6HTHETABE1 5 .8,9X1HOE15.8/4X6HODE G 80 
4LDXE15.8,6X4HDWDXE15.8,2X8HDI0DX( 2) El 5. 8, 2X8HDI 00X( 3 ) £15. 8 , 2X8HDI 0 G 81 
5DX(4)E15. 8/2X8HDI ODX ( 5 ) E15 . 8, 1X9HDSMALR8DXE 15 . 8/ ) G 82 

13 FORMAT ( 1H06X6HPRAT IOE 15* 8 , 2X6HQRATI 0E 15.8/) G 83 

END G 84- 
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CHER 


An additional test on the DELTN (5 0 ) iteration is performed in CHER. After 
PROPIT returns to CHER; KEZ is tested. If the RHO computed in THEP is too small, 
KEZ is set to 1. When this is the case, DELL = DELTN and CHER then calls MAIN to 
reinitialize the program with a new DELTN. CHER is called by DERSUB. The flow 
chart for subprogram CHER is as follows: 
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C 


1 


SUBROUTINE CHER 

DOUBLE PRECISION VVAR,CCUVAR 


COMMON /MAI NP/ AA1 (61 , AC, AO* AOVERB, A1 , 8BI ( 6 ) , CAPHO, CAPH1* C APQ*CC , C 
1C 1 1 CCI » CCU VAR. (8 I » Cl It C I MAX, COSTS, COS* * CWMT B* DA F t ODELDX* DOER < 8 ) f DD1 
2t DELL, DELS tDELT A f DELTN , DELU , DIODXC 6) *DI1DWC6) *DGSRB * DPI DW * DR8DX, DR 
301DW,DTBDX,DUSDWtDUiOWtDVSDWtDVlDWtD^DXtE6LEl{7)t EELE2(7) ,EELT(3) , 
4EE1 C6 ) , ERRV AL (7 ) , EMUREF, EOl * EPSR, EPS2 1 * GGQ( 6 ) , GG1 < 6 ) * HONE * ICGSW, I E 
5RR* I GEO , I TEXT , I TIES T *KEKQNT * KETEST *KKKK*GPERA * P * P EEP* PI N* PG, PSTAG , 
6PItQ0tRfRBtRBXf RHOI N, RHON, RHOO, RHOl , R 1 1 SINTB , SI NSQW * S I NW, $ MALLB , SM 
7ALRB tSTAENTtS SPEC, SWMTBtTtTCGtTEAN, THETAS, TO, T1,UIN,U I NStUOtUStUl, 
8V* VO* VS, VI ,VVAR(8) ,W, X, XI 0< 6) , XI 1 ( 6) , XMAX, ZO* 21 

VQNE=CCUV AR ( 1 ) 

CALL PROPIT CR,UINStXIOtTD*POtRHOOtITTESTtRHOIN,STAENT,ZDUM,KEZ,EP 
IS 21 * AO , VONEt C APHO , UO, PI N, EMUREF ) 

IF ( KEZ • EQ* 1 1 GO TO 1 

RETURN 

DELL^DELTN 

CALL MAIN 

END 


N 

N 

N 

N 

N 

N 

N 

N 

N 

H 

N 

N 

N 

N 

N 

N 

N 

N 

N 

N 

N 

N 


1 

2 

3 

4 

5 

6 

7 

8 
9 

ID 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 - 
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PROPIT 


PROPIT is a subprogram called by CHER during each integration step to solve for 
the body properties through a Newton-Raphson iteration scheme. PROPIT has a calling 
sequence in which R, UINS, RHOIN, PIN, ITTEST, EPS21, and EMUREF are constants. 
The variables XIO, AO, and CAPHO are computed for the particular value of VONE, the 
independent variable X, when PROPIT is called. The body properties, PDUM, RHODUM, 
UGEN, TGEN, and ZDUM and the static enthalpy STAENT are computed by the subpro- 
gram. An indicator KEZ is set in THEP when RHO is too small. The flow chart for 
subprogram PROPIT is as follows: 
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SUBROUTINE PROPIT { R ,UI NS , X 1 0, TGEN, PDUM ,RH0DUM, ITTEST , RHOI N , STAENT M 1 

1,ZDUM,KEZ,EPS21,A0,V0NE,CAPH0,UGEN,PIN,EMUREF) M 2 

C M3 

DIMENSION XIO (6) 

C KTEST IS THE NUMBER OF ITERATIONS ALLOWED ON RHO M 5 

KTE ST=0 M 6 

KEZ=0 M 7 

CAPHO=XIO(5)/XIO(2) M 8 

XI02S=XI0( 2 )**2 M 9 

IF (VONE-.005) 1,1,3 M 10 

1 UGEN=AO*VONE M 11 

UGENS=UGEN**2 M 12 

STAENT=CAPH0-UGENS/2. M 13 

RH0DUM=XI0(2)/UGEN M 14 

RUGS=RHODUM*UGENS M 15 

PDUM=XIO( 3J-RUGS M 16 

M 17 

M 18 

RETURN M 19 

STAENT=CAPHO- . 5*X I02S/ ( RH0DUM**2) M 20 

CALL THEP (R, PDUM,RHODUM, TGEN»ZDUM»STAENT,EMUREF,PIN,RHOIN»UINS»KE M 21 

1Z> M 22 

IF (KEZ.EQ.l) RETURN M 23 

FRODUM=RHODUM— ( XI 0( 3 > -PDUM ) * ( RH0DUM**2 ) / XI02S M 24 

DELRHO=— • 1 M 25 

RHODMP=RHODUM+DELRHO M 26 

STAD=CAPHQ— »5*XI02S/(RH0DMP**2) M 27 

CALL THEP (R,PD,RHODMP,TD,ZD, STAD, EMUREF,PIN,RHOIN, UI NS,KEZ > M 28 

FPRHO= ( RHODMP— (XI0(3)-PD)*( RH0DMP**2) / XI02S-FR0DUM) /O ELRHO M 29 

FROFPR=FRODUM/FPRHO M 30 

AFRFPR=AB S( FROFPR ) M 31 

IF I AFRFPR-EPS21 ) 6,6,4 M 32 

KTEST=KTEST+1 M 33 

IF (KTEST— ITTEST) 5,7,7 M 34 

RHO DUM=RHODUM— FROFPR M 35 

GO TO 3 M 36 

UGEN=XIO( 2 ) /RHOOUM M 37 

UGENS=UGEN**2 M 38 

RUG S=RHODUM*UGENS M 39 

GO TO 2 M 40 

WRITE (6,8) M 41 

KEZ=1 M 42 

RETURN M 43 

M 44 

M 45 

FORMAT ( IH06X34HI TTERATION LIMIT REACHED IN PROPIT) M 46 

END M 47 
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COEF 


Subprogram COEF is called by DERSUB to calculate the coefficients of the gov- 
erning differential equations. The flow chart for subprogram COEF is as follows: 
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SUBROUTINE COEF J 1 

C J 2 

DOUBLE PRECISION VVAR*CCUVAR J 3 

C J 4 


COMMON /MAI HP/ AA1 ( 6 ) * AC, AO, AGV ERB* A1 * BB1 { 6 ) * CAPHO* CAPH1, C APQ f CC * C J 5 
1C1,CCI*CCUVAR(8),CI I*CIMAX* COSTB f COSW,CWMTB* DAF*DDELDX* DDER(8) ,DD1 J 6 
2* DELL* DELS * DELTA* DELTN* DELU* DIODX ( 6) * D I 1DW ( 6 1 , DOSRB, DPIOW , DRBDX* DR J 7 
3Q1DW*DTBDX*DUSDW*DU1DW*DVSDW*DV1DW,DWDX*EELE1 (7),EELE2<7) *EELT<3) * J 8 
4EEK6) * ERRV AL (7)* EMUREF , EO 1 * EPSR , EPS 2 1 , GGO ( 6 ) * GG1 { 6 ) , HONE * ICOSW*IE J 9 
5RR* I GEO * I TEXT t.ITTEST *KEKGNT * KETEST * KKKK * OPERA* P* PEEP* PI N, PO* PSTAG* J 10 
6P1* QD* R*RB*R8X* RHOIN* RHON* RHOO* RH01* R1*SINTB*SINSQW*SINW*SMALLB*SM J 11 
7ALRB*STAENT*SSPEC*SWMTB*T,TCG*TEAN,THETAB,T0* T1*UIN* U I NS* UO* US*U1 , J 12 


8V*V0*VS*Vl,VVAR<8} *W*X*X10(6) *XI1 ( 6) ,XMAX,Z0*Z1 J 13 

C J 14 

CAPR1-R1 J 15 

C J 16 

T PT 0D=3 .+2**DAF J 17 

AA1<2)=3.*(XIGC2)-XI1( 2) )+2.*DAF*(X 10 ( 2 )-X 1 1 ( 2 ) ) J 18 

AA1(3)=3.*(XI0(3I-XI1(3))+2.*DAF*(XI0<3)-XI1< 31 I J 19 

AA1 (4)~— TPT0D*XI1 (4 ) J 20 

AA1<5)=3.*<XI0<5)-XIH5)) + 2.*DAF*(XI3( 5)-XIlC 51) J 21 

DO 1 J=2* 5 J 22 

B81(J)=DELTA*TPT0D*DI1DW(J> J 23 

1 CONTINUE J 24 

CC1=DELTA*(3.+DAF) J 25 

DD1-1. J 26 

EGl=-( l.+CAPQ*OELTAl/CAPQ*SyMTB/CWMTB*DTBDX J 27 

GG1(5)=RH01*V1*CAPH1 J 28 

GGI(2)=RHG1*V1 J 29 

GG1 ( 3 ) =RH01*U1*V1 J 30 

GGO (4 ) -PO J 31 

GG1(4)=P1+RHG1*V1**2 J 32 

EEl<2>^3.*D0SRB*CXia{2)+XIl<2))*DRBDX-DELS*0PERA*<XI0(2)+2.*XIl<2) J 33 
l)*DTBDX+6.*(l.+QD)*{i.+DAF)*GGl(2)+0ELTA*TPT00*GGi<2l*DT3DX J 34 

EE1 (3) =3.*D0SRB*C RH00*U0**2 +RHO 1*U 1**2 ) *DRBDX-DELS*OPERA* ( RHOO*UQ* J 35 
1*2+2. *RH01*U1**2) *DT8DX+QD*<3.+2.*DAF)*XU(4)+6.*(1*+QD)*C l.+DAF)* J 36 
2X1 1 C 4) +2* *DELTA*TPTOD*X 1114) *DTBDX J 37 

EE1(4)«3* *DGSRB*X II (4)*DRBDX-2. *DELS*OP ERA*X 1 1 ( 4)*DTBDX-3.*(QQ+DAF J 38 
ll*<XIG(3)+XIim)~2.*QD*DAF*(XIO{ 3)+2.*XU(3) )+6.*( 1 • +QD) *( 1 . +DAF ) J 39 
2*GG1 ( 4 )—6 • *GG0C 4 ) +QD*D AF*< RHOO*UO**2 + 2 **RHQ 1*U1**2 ) +3 ♦*DAF*( RHOO*U J 40 
30**2+RHOI*U1**2)+DELTA*{3.+2.*DAF) *RH01*( V1**2~U1**2) *DTBD X J 41 

6E1(5)=3**0QSRB*< XIO( 5>+XIl(5) ) *DRBDX-DELS*GPERA* < XI 0 ( 5 H-2 .*XI 1 ( 5 ) J 42 
l)*0TB0X+6.*(i.+QD)*U.+0AF)*GGl(5 ) +6 . *R1+DELTA* { 3 .+2. *DAF > *GG1 ( 5 ) * J 43 
2DTBDX J 44 

RETURN J 45 

END J 46- 


51 



CALD 


The derivatives of the governing differential equations which are used in DERSUB 
are computed in CALD. CALD is called by DERSUB. The flow chart for subprogram 
CALD is as follows: 
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SUBROUTINE CALD 0 1 

DOUBLE PRECISION VVARtCCUVAR 0 2 

D 3 
0 4 

COMMON /MAINP/ AA 1 (6 ) , AC, AO, AQVER8 , A1 , BB1 ( 6 > , CAPHO, CAPH1, CAPQ ,CC ,C D 5 
1C1,CCI ,CCUVAR<8> ,CI I,CIMAX,COST8,COSlrf,CWMTB,DAF,DDELDX,ODER{8) ,DD1 0 6 

2, DELL, DELS, DELTA, OELTN,OELU, DIODXI 6) ,DI 1DWI6) , DOSRB, DPI DW, DRBDX , DR D 7 
301DW,DTBDX»DUSDW» DU1DW, DVSDW, DV1DW* Di»lDX,EELEl(7)*EELE2(7),E£LT{3) , 0 8 

4EEK6) , ERRV AL ( 7 > , EMUREF ,£Q1,EPSR,EPS21» GGO( 6 ) » GG1 ( 6 ) , HONE , ICOSW, I E D 9 
5RR, I GEO, ITEXT, ITT EST,KEKONT, KETEST,K<KK, OPERA, P, PEEP, PIN, PO.PSTAG, D 10 
6P1,QD,R»RB,RBX» RHOI N, RHON, RHOO, RHOl , R1 , SINTB, SI NSQW , S I NW, S MALLB , SM D 11 
7ALRB*STAENT»SSPEC»SWMTB*T*TCG»TEAN,THETAB,T0,T1,UIN,UINS,U0*US»U1, D 12 


8V, VO, VS, VI, WAR (8) ,V),X,XIO(6) , X I H 6 ) ,XMAX,Z0,Z1 D 13 

0 14 

RBS=SMALRB**2 ' D 15 

RACE=RBS/(SMALLB**2-RBS) 0 16 

IF ( IGEO— 2 ) 3,1,2 D 17 

DRBDX=1 . /SQRT ( 1 .+ AO VERB** 2 *RACE ) D 18 

GO TO 3 D 19 

ORBDX=l»/{ 1 «+AOVERB**2*SMALRB/ { SMALL8**2+RBS>) D 20 

DDELDX=E01 D 21 

DWDX= ( - AA 1 { 4) *DDE LDX— EE 1(4) )/BBl(4) D 22 

IF ( DWDX ) 5,5,4 D 23 

DELU=DELTN D 24 

CALL MAIN D 25 

DI0DX(2)=(-AA1(2)*DDELDX-BB1(2)*DWDX-EE1<2})/CC1 D 26 

DIODX ( 3 )=UO*D IODX 1 2 I 0 27 

DIODX(5)=( -AA 115) *DDELDX-BB1 (5)*DWDX-EE1(5) )/CCl 0 28 

RETURN D 29 

END 0 30- 
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DERSUB 


DERSUB is a subprogram used by INTI, the integration routine, to evaluate the 
derivatives. The flow chart for subprogram DERSUB is as follows: 
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SUBROUTINE DERSUB H 1 

DOUBLE PRECISION VVAR,CCUVAR H 2 

C H 3 

COMMON /MAI NP/ AAI { 6) , AC, AO, AOVERB , AI , BBI ( 6 ) , CAPHO, C APH l, C APQ,CC , C H 4 
1C1,CCI»CCUVARI81,CI I ,C I MAX , COSTB, COSH , CWMT B, DAF, DDELD X, DDERI 81 ,DD1 H 5 
2, DELL, DELS, DELTA, DELTN,DELU,OIODXl 61 ,DI1DWI 61 ,OOSRB,OP1DW,DRBDX,OR H 6 
301DW,DTBDX,DUSDW,DUlDW,DVSDW,DVlDW,DtfDX,EELElI71,EELE2I 7) ,EELTI31 , H 7 
4EE1I61,ERRVALI71» EMUREF »ED 1 , EPSR, EPS2 1 » GGO( 6 1 ,GGI ( 6 1 , HONE ,ICOSW,IE H 8 
5RR, I GEO, I TEXT , ITTEST,KEKONT ,KETEST,X<KK, OPERA,?, PEEP, PIN,PO,PSTAG, H 9 
6P1»QD,R»RB»RBX,RH0IN» RHON, RHQ0,RH01,R1,SINTB,SINSQW,SINW»SMALLB,SM H 10 
7ALRB, STAENT,SSPEC , SWMT B,T,TCG,TEAN,THETAB,T0,T1,UIN,UINS, UO, US ,U1 , H 11 
8V,V0,VS,V1,VVAR(81,W,X,XI0(61,XI1(61,XMAX,Z0,Z1 H 12 

C H 13 

X=CCUVAR( 1 1 H 14 

OELTA=CCUVAR( 2) H 15 

W=CCUV AR( 3 1 H 16 

DO 1 1=2,5 H 17 

X 10 ( 1 1 = CCU VAR (1+21 H 18 

1 CONTINUE H 19 

IF ( IGE0.GE.2 1 SMALRB=CCUV AR ( 8 1 H 20 

CALL GEO H 21 

CALL RANH H 22 

DOSRB=DELT A/SMALRB H 23 

COSTB=COS (THE TAB 1 H 24 

SINTB=SIN(THETAB1 H 25 

DAF=DELTA*COSTB/SMALRB H 26 

OPERA=S INTB/SMALRB H 27 

CALL CHER H 28 

CALL COEF H 29 

CALL CALD H 30 

ODER! 2 1 =DOELO X H 31 

DDERI 3 1 =DWDX H 32 

DO 2 1=2,5 H 33 

DDERI I+21 = DlODXm H 34 

2 CONTINUE H 35 

IF IIGE0.GE.21 DDER(8)=0RBDX H 36 

RETURN H 37 

END H 38- 
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CHSUB 


CHSUB is called by INTI. It is written by the user and is described in appendix A. 
In program D1250 the features of CHSUB are not needed; therefore, CHSUB is a dummy 
routine. 


SUBROUTINE CHSUB 

RETURN 

END 


E 1 
E 2 
E 3- 


56 









o o 


FUNCTION FOFX IPQUM) 0 1 

0 2 
0 3 

DOUBLE PRECISION VVAR,CCUVAR 0 4 

COMMON /MAINP/ AA 1 ( 6 ), AC, AO ? AOV ERB , A1 , B BI ( 6 > , CAPHO , C A PH 1 , CAPQ , CC , C 0 5 

ICIfCCI fCCUVAR (8) ,CI I , C I MAX , CGSTB, CQSW , CWMT B, OAF , DDELDX, ODE R ( 8 ) ,DQ1 0 6 

2, DELL, DELS, DELTA, DELTN , DELU , DIOOX { 6 1 , D I 1DW(6J , DOSR8 , DP 1 DW , DRBDX, DR 0 7 

301DW,DTBDX, DU SOW, DU IDW , DVS DW, DV IDW , DW DX , EEL E 1 < 7 ) , EEL E 2 ( 71 , E E L T < 3 ) , 0 8 

4EE1I6) , ERRVALC7I, EMUREF ,EO 1 , EPSR, EPS2 1 , GGG ( 6 ) ,GG1(6) ,HONE, ICGSW,IE 0 9 

5RR, IGEO, I TEXT, ITT EST,KEKOMT, KETEST,K<KK, OPERA, P, PEEP, PIN, PO,PSTAG, 0 10 

6PI,QD,R,RB,RBX f RHOIN,RHON,RHOO,RHOI,Rl ,SINTB, S l NSQW , S I NW, S MALLB , SM 0 II 
7ALRB,STAENT ,SSPEC, SWMTB,T ,TCG,TEAN,THETAB,T0,T1,UIN,UINS, UO , US ,U1 , 0 12 


8V, VO, VS, VI, WAR (8) , W» X , XI Q ( 6 ) , X II ( 6 > , X M AX , ZO, Z I 0 13 

ROPE=- ( AC**2*(P1-PDUM)**2)/ ( PEEP*CC** 2*PDUM ) 0 14 

1 STAPE=CAPHl-3.*Rl/< AC*(Pl-PDUMn 0 15 

IF (STAPE.GT.O) GO TO 2 0 16 

Rl=*9*Rl 0 17 

GO TO 1 0 18 

2 CALL THEP I R ,POX, ROP£,TD, ZD, STAPE, EMUR EF, P IN, RHOIN, UI NS ,KEZ ) 0 19 

FOF X=PG X 0 20 

RETURN 0 21 

END 0 22- 
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PROD 


Subprogram PROD reads and writes the radiation NAMELIST input RADI and calls 
the subprogram to compute radiation. The thermodynamic properties within the shock 
layer are computed in the 77 distribution. PROD has a calling sequence in which USTAR 
determines whether the properties are off the stagnation point. If they are off the stagna- 
tion point, USTAR becomes the velocity at the body surface. IPROD is the option to read 
the radiation input data. AUODUM and AU1DUM represent ap and aj when PROD is 
called by MAIN and uq and uj when PROD is called by CONT. QRYO and QRY1 are 
computed by RAD FLUX and are the total radiation heat fluxes at the body and shock, 
respectively. The flow chart for subprogram PROD is as follows: 
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o non 


SUBROUTINE PROD ( USTAR , IPROD, AUODUM, AU1DUM, QRYO , ORY 1 ) P 1 

C P 2 

DOUBLE PRECISION VVAR,CCUVAR P 3 

C P 4 

COMMON /MAINP/ AAK6) , AC, AO , AOVERB, A1 , BB1 ( 6 ) , CAPHO, CAPH1, C APQ, CC, C P 5 
1C1,CCI»CCUVAR(8)»CII»CI MAX»COSTB» COS W » CWMT B» DAF *DDELDX»DDER(8) , DD1 P 6 
2,DELL,DELS,DELTA,DELTN,DELU,DIODX(6),DI1DW(6),OOSRB,DP1DW,DRBDX,DR P 7 
301DW»DTBDX»0USDW*DU1DW*DV5DW,DV1DW,DW0X»EELE1(7)»EELE2(7),EELT(3) , P 8 
4EE1 { 6 ) , ERRVAL (71, EMUREF , EDI , EPSR, EPS2 1, GGO( 6 ) , GG1( 6 ) , HONE , ICOSW, I E P <3 
5RR, IGEO*ITEXT»ITTEST ,KEKQNT, KETEST , KKKK, OPERA, P, PEEP, PIN, PO,P ST AG, P 1C 
6P 1 , QD, R ,RB, RBX, RHOI N, RHON, RHOO, RH01»R1»SINTB» SI NS3W, S I NW, SMALLB , SM P 1] 
7ALR8,STAENT,SSPEC» SWMTB ,T,TCG»TEAN»THETAB,TO,Tl,UIN»UINS,UO,US»Ul, P V, 
8V,VO,VS,V1,VVAR(8),W,X,XIO(6>,XI1(6),XMAX»ZO,Z1 P l: 

COMMON /WON/ QRYPC (20) , QRYPL( 20 ) P 1< 

COMMON /AAA/ A( 6, 6) , A12 ( 16 , 1 2 ) , A22 ( 7, 7 ) , ALP ( 6 ) , ALPT< 1 6 , 5) , BMT< 16) , P 1' 
1C (6,16) ,CH(16,2),CP(16) ,DEGIC6),H( 16 ) , J AT ( 16 , 5 ) , JPH < 1 6 ) , KA T( 6 ) ,KOD P If 
2E ( 1 6 ) , RA( 16,2) ,RBO( 16 , 2 ) ,RC ( 16 , 2 ) , RDC 16 ,2 ) , RD1 ( 16 ,2 ) , RE ( 1 6 , 2 ) , RE1 ( P 1’ 
316,2), SD( 16), TB(3)»TU( 16,2 ) ,TU2( 16,2) ,VN( 17) , VNE( 16) , VNT(17),VNU( 1 P li 
46,6) ,W26( 6) » W3( 16) ,Y( 16) »RF (16,2) » RC1 ( 16,2 ) ,HS , IG , I GMS , I GMS P , IGP, I P 1< 
50N, IS , I SP , I SPNGP, ISPNG2 ,N, NG »NP, PRESS, W27,R HO, WM,SYU, I SPNG P 2C 

COMMON /RAD/ YY( 100 ) , TE E( 1 00 ) , FHVt 20 ) , NHV, NY , C2 , I Y P 21 

COMMON /RAD/ XNN(7,100) P 22 

COMMON /RAD/ NIHVC , FHVC ( 50 ) , AHV (50) , AH VL( 20 ) P 23 

COMMON /RAD/ Cl ,C3, C4, FLG, C5 , FLG1 P 24 

COMMON /RAD/ YDELT, DELTD, FL1 , FL2 P 25 

COMMON /RAO/ GEE( 8 > , EP S ( 8 ) , NU( 20 ) , HOI 70 ) , FF ( 70) ,GAMP ( 70 ) , WOL ( 20) , F P 26 
1HVMI20) »FHVP( 20) P 27 

COMMON /CIONCL/ F C 100 , 10) , F2 < 100, 10 ) , HVL( 70 ) , EP ,K2, K1 , I FL , I YCON, I Q P 28 
11 *WMI ,BIJ( 100,10) ,GMIN( 100 ,10), GPLU (100*10), IAED P 29 

DIMENSION NICN(IOO) P 30 

DIMENSION HH(IOO), PRES(IOO) P 31 

REAL KN*KNY P 32 

NAMELIST /RADI/ GEE, EP S ,FHVM, FHVP , FHV , WOL , NU , ND ,HVL , FF , GAMP, I S, N I H P 33 
1VC,FHVC»AHV»AHVL»FL1»FL2,FLG»FLG1,NIC»NICN» IOPT,TB ,NY P 34 

P 35 

THIS ROUTINE IS CALLED IN THE MAIN PROGRAM WHEN IPR0D=1 TO P 36 

READ IN THE DATA FOR RADFLUX,WHEN IPROD =2 , NAMELIST IS NOT CALLED P 37 

P 38 

IF ( IPROD— 1 ) 2,1,2 P 39 

1 READ (5, RADI) P 40 

WRITE (6, RADI) P 41 

WRITE (6,8) P 42 


60 



vO oo n o ^ o ui ^ 


2 UINS=UIN**2 P 43 

I Q=NY P 44 

YY(1)=0 P 45 

IF IUSTAR) 3,3,4 P 46 

3 HH(1)=CAPH0*UINS/4,19E7 P 47 

GO TO 5 P 48 

HH ( 1) = ( CAPHQ-UQ**2/2 . )*UINS/4. 19E7 P 49 

PRES ( 1 )=P0*RH0IN*UINS/ 1.01325E6 P 50 

IF (KEKONT .EQ.KETEST) WRITE (6,9) P 51 

ETA=0 P 52 

STAENT =CAPH0— 110**2/ 2 . P 53 

IF (KEKONT. EQ.KETEST) WRITE (6,10) ETA , PO, STAENT , CAPHO, AUODUM, VO P 54 

DO 6 N=2 , NY P 55 

KN=N P 56 

KNY=NY P 57 

ETA=( KN— 1 • ) /( KNY-1 • ) P 58 

YY(N) = ETA P 59 

RHO AU= RHOO* AUODUM +( RHO 1 *AU I OUM-RHOO* AUODUM ) *E T A P 60 

RHOAUS=RHCO*AUODUM**2+(RHQ1*AU1DUM**2-RHQO*AUODUM**2) *ETA P 61 

RHQAUH=RHOO*AUODUM*CAPHO+( RHOi*AUiDUM*CAPHl-RHOO*AUODUM*CAPHO )*ETA P 62 

RH0AUV=RH01*AU1DUM*V1*ETA P 63 

AUDUM=RHOAUS/RHOAU P 64 

V=RHOAUV/RHOAU P 65 

IF (USTAR.GT.O) UST AR=AUDUM P 66 

SUMH=RHOAUH/RHOAU-USTAR**2/2.-V**2/2. P 67 

STAENT® SUMH P 68 

CAPH=RHOAUH/RHOAU P 69 

PDUM=PO+(QD*RHOO*UO**2)*ETA+(P1-P3-QD*RHOO*UO**2)*ETA**2 P 70 

IF (KEKONT. EQ.KETEST) WRITE (6,10) ETA, PDUM, STAENT, CAPH, AUDUM, V P 71 

HH( N)=SUMH*UINS/4 • 19E7 P 72 

PRE S( N )=PDUM*RHOI N*UI NS/1, 01325E6 P 73 

CONTINUE P 74 

IF (KKKK.EQ.O) GO TO 7 P 75 

CALL RADFLUX ( N IC ,N ICN , IOP T , PRES, HH, QRYO, QRY1 ) P 76 

IF (KEKONT. EQ.KETEST) WRITE (6,11) QRYPC( 1 ) » QRYPC ( N Y ) »QRYPL( 1) ,QRY P 77 

1PL( NY) ,QRY0,QRY1 P 78 

RETURN P 79 

P 80 
P 81 

FORMAT < 1 HI I P 82 

FORMAT ( 1H04X31HD I STRI BUTION ACROSS SHOCK LAYER/12X3H ETA14X4HPDUM1 P 83 

13X6HSTAENT1 1X4HCAPH13X5HAUDUM12X1HV/J P 84 

10 FORMAT (6 ( 2XE 15.8 ) ) P 85 

11 FORMAT ( 1H0 ,6X» 5HQRYPC2E15 . 8»2X»5HQRYPL2E15.8»2X» 4HQR YOE 15 . 8, 2X » 4H P 86 

1QRY1E15 .8) P 87 

END P 88- 
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RADIATION SUBPROGRAMS 


The radiation package was adapted from program RATRAP (ref. 4). Reference 3 
describes RATRAP as it is applied to D1250. If the input quantity KKKK is greater than 
zero, a radiation solution is computed. Subprogram PROD links D1250 to RATRAP by 
reading the input for the radiation program and by calling RAD FLUX. RAD FLUX serves 
as the main program for the radiation subprograms. A listing of the radiation subpro- 
grams follows. 



o n 


c 


SUBROUTINE RADFLUX C NIC, NICN, IOPT, PRES , HH, QRYO, QRY1 ) Q 

Q 1 

DOUBLE PRECISION VVAR.CCUVAR Q 3 

Q 4 
Q 5 

COMMON /MAINP/ AAU6) , AC, AO, AOVERB, A1 , BBI < 6 ) , CAPHO, CAPH1, C APQ, CC , C Q 6 
1C1,CCI,CCUVARC8),CII,CIMAX,C0STB,C0SW,CWMTB,DAF,DDELDX,DDER(8),DD1 Q 7 
2,DELL»D£LS»DELTA, DELTN»DELU»DI0DXC6), DI 1DW ( 6 ) » DOSRB, DP 1DW* DRBDX, OR 3 8 

301DW,0TBDX,DUSDW,DU1DW,DVSDW,DV1DW,DWDX,EELE1C7),£ELE2C7),EELT(3> , Q 9 
4EE1C6) , ERRVALC7), EMUREF , EDI , EPSR, EPS2 1 , GGO t 6 > , GG1C 6 ) , HONE , ICOSW, I E Q 10 
5RR, IGEO,ITEXT,ITTEST,KEKONT,KETEST,KKKK,OPERA,P,PEEP, PIN, PO,PSTAG, Q 11 
6P 1 » QDt R*R B,RBX, RHOIN, RHON, RHOO, RHO 1, R 1 , SINTB, SINSQW,SINW,SMALLB,SM Q 12 
7ALRB,STAENT,SSPEC,SWMTB,T,TCG,TEAN,THETAB,T0,T1,UIN,UINS,U0,US,U1, Q 13 
8V,VO,VS,V1,VVARC8>,W,X,XIOC6),XI1<6),XMAX,ZO,Z1 Q 14 

COMMON /WON/ QRYPC (20) ,QRYPL(20) Q 15 

COMMON /AAA/ AC 6, 6) , A12 C16 , 12 ) , A22 ( 7, 7 ) , ALP( 6 ) , ALPT ( 16 , 5) , BMT ( 16) , Q 16 
1C (6, 16) ,CH(16»2)»CP(16),DEGI(6)»H{ 16 ) , JAT ( 16, 5 ) , JPH ( 16 ) ,KAT< 6 ) ,KOD Q 17 
2EC16) ,RA( 16,2),RB0(16,2),RC( 16,2) ,RD< 1 6 ,2 1 , RD1C 16 , 2 ) , RE ( 16 , 2 ) , RE1 ( Q 18 
316,2),S0( 16),TB(3),TU< 16,2) ,TU2C16,2) ,VN(17) , VNEC 15 ) , VNTC 17 ) , VNUC 1 Q 19 
46,6),W26(6),W3(16),Y(16),RF(16,2), RC1C 16,2) ,HS, IG, I GMS , IGMSP , IGP, I Q 20 


50N, IS, ISP,ISPNGP,ISPNG2,N,NG,NP,PRESS,W27,RH0,WM,SYU, ISPNG Q 21 

COMMON /RAD/ YY < 100 ) , TEEC 100 ) , FHVC 20 » , NHV, NY, C2 ,1 Y Q 22 

COMMON /RAD/ XNNC7,100) Q 23 

COMMON /RAD/ NIHVC, FHVC C 50 ) , AHV C 50 > , AH VLC 20 ) Q 24 

COMMON /RAD/ Cl ,C3, C4, FLG, C5, FLG1 Q 25 

COMMON /RAD/ YOELT, OELTD, FL 1 , FL2 Q 26 

COMMON /RAD/ GEEC 8 ) , EPS ( 8 ) , NUC 20) , ND( 70 ) ,FF ( 70) ,GAMP ( 70 1 , WOL( 20) ,F 3 27 

1HVMC20) ,FHVP( 20) Q 28 

COMMON /CIONCL/ F ( 100, 10) , F2 < 100, 1 0 ) , H VLC 70 > , EP ,K2, K1 , I FL, I YCON, I Q Q 29 

ll,WMI,BIJ(100fl0) ,GMIN( 100 , 10) , GPLUC 100 ,10 ) , I AED Q 30 

DIMENSION NICN(IOO) 0 31 

DIMENSION PRES(IOO) Q 32 

DIMENSION HHC100) Q 33 

NHV=1 8 Q 34 

C 4=1. 273 Q 35 

G=1 . Q 36 

N ICNC 2 ) =NY Q 37 

IF (I0PT.EQ.2) GO TO 1 Q 38 

TBC 1 )=TEE ( 1 ) Q 39 

DO 15 1=1, NY Q 40 

PRE SS=PR£ SCI) Q 41 

IF CI0PT.EQ.2) 60 TO 2 Q 42 
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TB(2)=TB<1) Q 43 

GO TO 3 Q 44 

2 TB(2) = TB( 1)*1 .2 Q 45 

HS=HH(I» Q 46 

3 I 1=0 0 47 

IF (I.EQ.1J GO TO 5 Q 48 

4 Z=29.*PRESS/( 82.057*TBI2J*RHO) Q 49 

IF (Z.LE.1.5) GO TO 6 Q 50 

5 CALL HTP (IZW) Q 51 

IZW=IZW Q 52 

GO TO 7 0 53 

6 CALL LTP (IZW) Q 54 

IZW=IZW Q 55 

7 CALL FEMP (IHELP.IOPT) g 56 

IF (I0PT.EQ.2 I 60 TO 8 Q 57 

TB( 2)=TEE ( I ) g 58 

TB( 1)=TB( 2) Q 59 

GO TO 9 g 60 

8 TEE ( I >=TB ( 2 ) Q 61 

HS=HH ( I ) g 62 

TB( 1 )=TB( 21/1.2 g 63 

9 IF (IHELP.EQ.O) GO TO 10 Q 64 

WRITE ( 6t 17 ) I»PRESSyHS,TB( l)yTB(2) Q 65 

STOP g 66 

10 IF (Il.NE.O) GO TO 11 Q 67 

11=1 g 68 

GO TO 4 g 69 

11 IF I IZW. EQ. 21 GO TO 12 g 70 

XNN ( 1 » I )=0 g 71 

XNN ( 2y I > = FLUT ( BMT t.3 ) • W3 ( 3 1 » RHO ) 0 72 

XNN (3*1 )=FLUT(8MT{2)»W3{2)t RHO) Q 73 

XNN I 4» 1 )=0 • Q 74 

XNN I 5» I )=0. Q 75 

XNN ( 6 » I )»FLUT(BMT,W3,RH0J g 76 

XNN (7*1 )=FLUT (BMT (4) ,W3(4) ,RHO) Q 77 

GO TO 14 g 78 

12 DO 13 L = 1 » 4 Q 79 

13 XNN ( Ly I ) = FLUT I BMT tL+2 ) » W3( L+2 ) » RHO ) Q 80 

XNN ( 5y I )*FLUT(BMT,W3yRH0) g 81 

XNN (6y I )=0. Q 82 

XNN (7, I )=FLUT ( BMT (2) ,W3I2) ,RH0) g 83 

14 TB( 2)*T EE ( I ) g 84 

RH0ND=RH0/RH0 IN g 85 


64 



o o 


IF (KETEST.NE. KEKGNT ) GO TO 15 Q 86 

IF (M0D(I-1,50).EQ.0) WRITE (6,18) Q 87 

WRITE (6,19) YY (I ),PRES(I),TEE( I) , HS, RHO, RHONO Q 88 

15 CONTINUE Q 89 

C2=OELTD*FLl Q 90 

C1=DELTD*FL2 Q 91 

C3=C1/3.1416 Q 92 

C5=C1 Q 93 

EP=.001 Q 94 

DO 16 L=1,NIC 0 95 

IY=NICN(L) Q 96 

I YCON= I Y q 97 

YDELT=YY ( I Y ) Q 98 

CALL CONTM Q 99 

CALL LINE g ioO 

16 CONTINUE Q 101 

QRYO=-( 1 » E7*( QRYPC ( 1 ) +QRYPL ( 1 ) ) /( RH0IN*UIN**3 ) ) Q 102 

QRY1=1.E7*(QRYPC(NYH-QRYPL(NY) )/(RHOIN*UIN**3) 9 103 

RETURN Q 104 

Q 105 
Q 106 

17 FORMAT ( 10H1FEMP BLEW, 110, 4E20.6) g 107 

18 FORMAT ( 14H0 PATH LENGTH, 4X, 8HPRESSt)RE,3X, 11HTEMPERATURE , 4X, 8 HENT Q 108 

1HALPY»5X, 7HDENSITY, 7X, 1 1HD£NSITY( ND),/,4X9H(Y/DELTA),9X,5H(ATM),7X Q 109 
2,7H(0£G K), 6 X, 8 H{ CAL/GM ),5X,7H(GM/CC) ) g HO 

19 FORMAT (6X6E13.5) Q m 

END g H 2 - 
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SUBROUTINE LINE R 1 

C R 2 

C R 3 

COMMON /WON/ QRYPCI20) ,QRYPL<20) R 4 

COMMON /AAA/ A< 6, 6 ) , A12 ( 16 , 1 2 ) , A22 ( 7, 7 ) , ALP( 6 ) , ALPT ( 16 , 5) , BMT( 16 ) » R 5 
1C ( 6 , 1 6 ) ,CH(16,2),CP(16) ,DEGI<6) ,H< 16 ) , J AT ( 16 , 5 ) , JPH ( 16 > , KAT { 6 ) ,KOD R 6 
2E(16),RA( 16,2),RB0<16, 2),RC( 16,2),RD< 16,2) , RD1< 16,2 ) , RE < 16, 2 ) ,RE1 < R 7 
316,2) , SD( 16) , TB(3) ,TU( 16,2 } , TU2 ( 16,2 ) t VNC 17 ) , VNE 1 16 ) » VNTI 17) , VNU( 1 R 8 
46,6),M26(6),W3(16),Y(16),RF(16,2),RCi( 16,2) ,HS, IG, I GM S , I GMSP , IGP,I R 9 
50N, IS, I SP ,1 SPNGP, ISPNG2 ,N, NG, NP, PRESS, W27,RHO,WM,SYU, ISPNG R 10 

COMMON /RAD/ YY tlOO ) ,TEE( 100 ) , FHV { 20 ) , NHV, NY ,C2 , I Y R 11 

COMMON /RAO/ XXN(7,100) R 12 

COMMON /RAD/ NOD, FDO( 50 ), AHV ( 50 ) , AHVL ( 20 ) R 13 

COMMON /RAD/ Cl ,C3 , C4 , FLG, C5 » FLG1 R 14 

COMMON /RAD/ YOELT ,0ELTD»FL1»FL2 R 15 

COMMON /RAD/ GEEI 8 ) , EPS (8 ) , NU ( 20 ) , ND( 70 ) » FF ( 70 ) ,GAMP ( 70 ) , WOL< 20) , F R 16 
1HVM(20) ,FHVP(20) R 17 

COMMON /CIONCL/ F ( 100 , 10) , F2 ( 100, 10 ) , HVL< 70 ) , EP,K2, K1 , I FL , I YCON, I Q R 18 
11 »WMI »BI J( 100,10) ,GMIN( 100,10) ,GPLU( 100,10 ) , IAEO R 19 

DIMENSION NICN(IOO) R 20 

DIMENSION PRES ( 100) R 21 

DIMENSION HH( 100) R 22 

DIMENSION BEE(IOO), WWM(IOO), WPI100), FMU(IOO), T2I100), SS(IOO), R 23 
1 TAU(IOO) R 24 

DIMENSION FX1I20), FY1(20) R 25 

DIMENSION WA(IOO), WS(IOO) R 26 

DIMENSION FX( 20) , FY(20) R 27 

DIMENSION IWP(IOO), IWM(IOO) R 28 

I YCON= I Y R 29 

G= 1 . R 30 

C9=2.42E14 R 31 

C 11=1. /2. 42E14 R 32 

FIP=0. R 33 

FIM=0. R 34 

K2=0 R 35 

00 21 K= 1 , NHV R 36 

FX(K)=0. R 37 

FY(K)=0. R 38 

FX1 (K)=0. R 39 

FY1(K)=0. R 40 

K1=K2+1 R 41 

K2=K2+NU( K) R 42 
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DO 1 1=1, NY 

R 

43 

Tl=TEEl I)*8.62E-5 

R 

44 

BEE ( I ) =0 • 

R 

45 

QH=FHV(K»/T1 

R 

46 

IF (QW.LT.85.1 BEE! I )=5040 • *FHV ( K 1 **3/ ( EXP ( QW )-l. ) 

R 

47 

WWM C 1 1 = 0, 

R 

48 

HP ( I )=0 • 

R 

49 

CALL MU (FHV(K),T1, XXN( 1,1) »FMU( 1)1 

R 

50 

CALL TRAP (NY, YY,FMU,C2 ,TAU > 

R 

51 

X1=FHVM( K ) 

R 

52 

X2=FHVP(K) 

R 

53 

QF=FLOAT ( NU (K ) ) 

R 

54 

DO 16 J=K1,K2 

R 

55 

I J= J-Kl+1 

R 

56 

T2( 11=0 • 

R 

57 

ssm=o. 

R 

58 

F2C 1,1 J)=0. 

R 

59 

B I J ( 1 , I J ) =0 *0 

R 

60 

DO 5 1=1, NY 

R 

61 

T1=TEE( I )*8.62E— 5 

R 

62 

E1=1.-EXP(-FHV(K)/T1) 

R 

63 

J1=ND( J J 

R 

64 

IF (J1.EQ.7) GO TO 2 

R 

65 

PN=XXN(2, I)+XXN(4,I) 

R 

66 

FN=PN/(4.+10.*EXP 1— 2.384/T1 J+6.*EXP(-3.576/Tl) ) 

R 

67 

GO TO 3 

R 

68 

PNP=XXN(1»I l+XXNC 3,1) 

R 

69 

FN=PNP/ (7„+5,*EXP(— 1#9/T1) I 

R 

70 

FN=FN*GEE(J1)*EXP(-EPS( J1»/T1J 

R 

71 

S1=FN*FF( J)*.0266*E1 

R 

72 

GAM={ GAMP ( J)*G*XXN( 5,1) *( Tl**.25) «-l . 0E-06) *C9 

R 

73 

IF (I.EQ.l) GO TO 4 

R 

74 

DY=„5*( YY(IJ-YYCI-l)) 

R 

75 

T2( U = T2( I— 1 ) +DY* (S1/GAM+0S1/0GAM) *C3 

R 

76 

SS( I) = SS( I— 1 ) +0Y*( Sl+OS 1)*C 1 

R 

77 

F2( I, I J l=F2( I— 1, I J I+DY*(S1*GAM+0S1*0GAM)*C5 

R 

78 

BI J (I , I J) =BI J ( 1-1 , I JI+DY*( GAM+OGAMI *4 . 135E-1 5 

R 

79 

0S1=S1 

R 

80 

OGAM=GAM 

R 

81 

CONTINUE 

R 

82 

IF (IY.EQ.l) GO TO 11 

R 

83 

DO 8 1=1, IY 

R 

84 

FINT=T2(IY1-T2( I) 

R 

85 
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IF (FINT.GT.C4I GO TO 6 

R 

86 


WSU) = (SS(IY)-SS( I ) )*Cil 

R 

87 


F ( I , I J )=0 • 

R 

88 


GO TO 7 

R 

89 

6 

WS( I l=2.*SQRT(F2(IY,IJ)-F2( I ,IJ>)*C11 

R 

90 


F(I,IJ)=1. 

R 

91 

7 

mwm( n=wwMt n+wsi t> 

R 

92 


IF (I.EQ.IY) GO TO 9 

R 

93 

8 

GMINIIt U) = ((BU( IY,IJ)-BIJ(I»IJ) ) / ( Y Y ( IY)-YY(I) J )**2 

R 

94 


GO TO 10 

R 

95 

9 

GMI N( I * I J )=0 

R 

96 

10 

IF ( IY. EQ.NY) GO TO 16 

R 

97 

11 

DO 15 1=1 Y»NY 

R 

98 


FINT=T2<I )— T2 1 IYI 

R 

99 


IF (FINT.GT.C4) GO TO 12 

R 

100 


WA ( 11= ( SSI I)-SS(IY) ) *C1 1 

R 

101 


F U , I J ) =0 . 

R 

102 


GO TO 13 

R 

103 

12 

WA( I )=2.*SQRT (F2( I, I J >-F2 ( IY, I J ) ) *C 1 1 

R 

104 


F ( I,IJ )=1. 

R 

105 

13 

WP{ II = WP{ I J+WAI I) 

R 

106 


IF (I.EQ.NY) GO TO 14 

R 

107 


GPLUII , IJ)=((BIJ( I, IJ)-BIJ(NY»IJ) }/(YY( I)-YY(NY) ) 1**2 

R 

108 

< 

GO TO 15 

R 

109 

14 

GPLUII^ IJ>=0 

R 

110 

15 

CONTINUE 

R 

111 

16 

CONTINUE 

R 

112 


IF (IY.EQ.1) GO TO 18 

R 

113 


I FL = 1 

R 

114 


DO 17 1=1 » I Y 

R 

115 


IWM ( I )=0 

R 

116 


IF (WWM(I)/QF.LT.WOL<Kn GO TO 17 

R 

117 


IWM ( I )=1 

R 

118 


101 = 1 

R 

119 


IAED=1 

R 

120 


CALL HAFACE ( 1,X1,X2> 

R 

121 


CALL HAFACE 12,0. ,0.) 

R 

122 


WWM ( I ) = WM I 

R 

123 

17 

WWM ( I )=WWMC I ) *EXP (TAU( I I-TAUI IY ) ) 

R 

124 


CALL TRAP 1 (IY,WWM,BEE,-1.,FX(KII 

R 

125 


FX1 (K)=WWMC1» 

R 

126 


FX(K)=FX(K)*AHVL(K) 

R 

127 


F IM=FI M+FXIK) 

R 

128 
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18 N1=NY— IY+1 R 129 

IF (Nl.EQ.l) GO TO 20 R 130 

I FL=-1 R 131 

DO 19 I=IY,NY R 132 

IWP(I)=0 R 133 

IF IWP( I)/QF.LT.WOL(K) ) GO TO 19 R 134 

IWP(I)=1 R 135 

IQ1=I R 136 

I AED=2 R 137 

CALL HAFACE (1,X1,X2» R 138 

CALL HAFACE (2,0.,0.) R 139 

WP(I)=WMI R 140 

19 WP( I)=WP< I)*EXP(TAU(IY)-TAU( I).) R 141 

CALL TRAP1 (Nl,VIPt IY) ,BEE( IY),1.,FY(K ) ) R 142 

FY1 (K) = WP(NY) R 143 

FY(K)=FY( K)*AHVL(K) R 144 

FI P=FI P+F Y( K) R 145 

20 CONTINUE R 146 

21 CONTINUE R 147 

QRYPLI IY)=FIM+FIP R 148 

RETURN R 149 

END R 150- 
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fSi CO 


SUBROUTINE CONTM S 1 

C S 2 

COMMON /WON/ QRYPCI20) »QRYPLI 20) S 3 

COMMON /AAA/ A { 6, 6 ) , A12 ( 16 , 12 1 , A22 ( 7, 7 ) , AL P( 6 ) , ALPT ( 1 6 , 5) , BMT ( 16 ) , S 4 
1CI6»16)»CHI16»2)»CPI16)»DEGII6)»H( 16) , J ATI 16,5) , JPH ( 1 6 ) ,KAT ( 6) ,KOD S 5 
2E(16),RAI 16,2),RB0116,2),RC116,2) ,RDl 16 ,2 ) , ROl 1 16, 2 ) , RE 1 16, 2 ) , RE1 1 S 6 
316, 2), SOI 16) ,TB(3), TUI 16,2), TU21 16,2) , VNI 17) ,VNEI 16) , VNTI 17) ,VNUI 1 S 7 
46,6),W26I 6),W3l 16) ,Yl 16 ) ,RF ( 16, 2) ,RC1 1 16,2 ) ,HS, IG, I GM S , IGMSP , IGP, I S 8 
50N, IS , I SP , I SPNGP, ISPNG2, N,NG,NP, PRESS, W27, RHO, WM, SYU, ISPNG S 9 

COMMON /RAD/ YYI 100 ) ,TEEI 100) , DUN 1 20) , NOV,NY,C2 , I Y S 10 

COMMON / RAO/ XXN(7,100) S 11 

COMMON /RAO/ NHV, FHVI 50 ) , AHVI 50 ) , AHVLI 20 ) S 12 

COMMON /RAD/ Cl ,C3, C4 , FLG, C5, FLG1 S 13 

COMMON /RAD/ YDELT ,DELTD»FL1,FL2 S 14 

COMMON /RAO/ GEE! 8 ) ,EP S 1 8 ) , NUI 20 ) , NO! 70) , FF I 70) , GAMP I 70 ) , WOLI 20) , F S 15 
1HVMI20) »FHVP( 20) S 16 

COMMON /CIONCL/ F 1 100 , 10) , F2I 100, 10 ) , HVL 1 70 ) , EP , K2, K1 , 1 FL , I YCON, I Q S 17 
11,WMI,8IJ( 100,10) *GMINl 100 ,10),GPLU(100,10),IAEO S 18 

DIMENSION NICNUOO) S 19 

DIMENSION PRESI 100) S 20 

DIMENSION HHUOO) S 21 

DIMENSION FMUI100), TAUUOO), FIMI40), FIPI40), BEE 1 100 ) , EMUOO), S 22 
1 EEPI100), FIIPI40), FIIMI40) S 23 

DO 9 K = 1 » NHV S 24 

DO 1 1=1, NY S 25 

T 1= TEE I I ) *8.6 2E-5 S 26 

CALL MU (FHVIK),T1,XXN( 1,1 ),FMUII ) ) S 27 

BEE I I ) = 0 . S 28 

AQ=FHVIK)/T1 S 29 

IF (AQ.LT.85.) BEE I I ) = 5040. *FHV I K) **3/ I EXP I AQ )-l. ) S 30 

1 CONTINUE S 31 

CALL TRAP INY,YY,FMU,C2,TAU) S 32 

IF l IY.NE.1) GO TO 2 S 33 

FIMIK)=0. S 34 

GO TO 4 S 35 

DO 3 1=1, IY S 36 

EMI I)=l. -EXPIT AUI I)— TAUIIY) ) S 37 

CALL TRAP1 1 1 Y, EM , BEE ,— 1 . , FI MI K ) ) S 38 

FIMIK)=FIMIK)*AHV(K) S 39 

4 N1=NY-IY+1 S 40 

IF IN1.NE.1) GO TO 5 S 41 

F IP I K) =0. S 42 
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m o 


GO TO 7 S 43 

DO 6 I=IY f NY S 44 

EEP(I)=1»— EXP(TAUIIY)-TAUII)) S 45 

CALL TRAP1 (N1»EEP( IYIt BEE ( IY)tl.»FIP(KII S 46 

FIP(K)=FIP(K)*AHV(K) S 47 

7 IF (K.NE.ll GO TO 8 S 48 

FlIM(l)=0. S 49 

F I IP ( 1 )=0 * S 50 

GO TO 9 S 51 

8 CALL TRAP IK ,FHV,FIM, 1 . ,FI IM» S 52 

CALL TRAP IK,FHV,FIP,1.,FIIPI S 53 

9 CONTINUE S 54 

QRYPC C I Y) =F 1 1 MINHV I+FI IP( NHV ) S 55 

RETURN S 56 

END S 57- 
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SUBROUTINE MU I HV » XKT, XM,XAPNU ) T I 

DIMENSION XKT(l), XNN(l), XNOI1), XNH1I, XN02I1), XNN2U), XAPNUI T 2 

II. 1), XM( 1 1 T 3 

SQA=7. 2 5E-16 T 4 

XNN 1 1 ) = XM (4 1 T 5 

XNOI 1 ) = XM (2 I T 6 

XNI(1) = XM<1HXM(3) T 7 

XN02 ( 1 ) =XM ( 6 1 T 8 

XNN2I 1 )=XM( 7) T 9 

L=1 T 10 

1=1 T 11 

XMQL=1 » T 12 

BQ7=0. T 13 

CALL ZHV (HV»ZO»ZN*ZI,ZC) T 14 

XN=14.3/XKT(I) T 15 

XO= 13.4/XKTII) T 16 

XI=25« 5/XKT ( I ) T 17 

XX=HV/XKT ( I ) T 18 

E Q1 =SQ A*XNN( I )*XKT ( I ) *4 «5*EXP (— XN + XX ) *ZN/HV**3 T 19 

EQ2=SQA*XNO< I)*XK Til) *0 .88888889*EXP< -XO+XX ) *ZO/HV**3 T 20 

EQ3=SQA*XNI(I)*XKTm*4.0*1.33*EXP(~XI+XX)*ZI/HV**3 T 21 

BQ5=XM0L*XN02 ( I )*400.0*SQRT ( TANH( 0 . 0975/XKT ( I ) ) ) *EXP (-TANH 1 0. 195/ ( T 22 

12.0*XKT{I)))*((HV-8.56)/0.805)**2)/2.687E+19 T 23 

BQ6=XM0L*XNN2 ( I) *1. 2E- 1 7*EXP ( -A8S( HV- 13.6+1 1. 0-0. 603/ XKT(I) ) 1**1. 3 T 24 

1) T 25 

EQ4=EQ1*EXP(4.22/XKTI I )— XX) T 26 

EQ5=EQ2*EXP(4.22/XKTm-XX) T 27 

EQ7=EQ3*EXP( 11.2/XKK I )-XX) T 28 

EQ9=EQ5+XN0( I )*3. 6E-17/ (9.0+5.0*EXP I -1 .98/XKTI I > )+EXP(-4. 18/XKT(Ii T 29 

1)) T 30 

EPC=I.O-EXP(-XX) T 31 

IF (HV-4.22I 1,1,2 T 32 

1 XAPNUIL,I )=EPC*( EQ1+EQ2+EQ3 ) T 33 

GO TO 17 T 34 

2 IF (HV-10.8) 3,3,8 T 35 

3 XAPNU(L,I )=EPC*( EQ4+EQ5+EQ3 ) T 36 

IF CHV-7.0) 17,4,4 T 37 

4 IF IHV-9.2) 5,5,6 T 38 

5 X APNUI L , I )= XAPNUI L, 1 1 +BQ5+BQ7 T 39 

GO TO 17 T 40 

6 IF (HV-10.0) 7,7,17 T 41 

7 XAPNUI L, I > = XAPNU(L,IH-BQ7 T 42 
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GO TO 17 T 43 

8 IF IHV-12 .0 ) 9,9,11 T 44 

9 XAPNUI L»I )=EPC* I EQ4+XNN I I ) *5. 16E-1 7*EXP l-XN+10. 8/XKTI II )/l 4.0+10.0 T 45 

1*EXP 1-2.38/XKT II)) +6.0*EXP 1-3.57/XKT II)) I+EQ5+EQ7) T 46 

IF IHV-11.0) 17,10,10 T 47 

10 XAPNUIL,! )=XAPNUIL,I)+BQ6 T 48 

GO TO 17 T 49 

11 IF (HV— 13.41 12,12,13 T 50 

12 XAPNUI L , I )=EPC*( EQ4+I 5 • 16E- 17*XNNl I ) *EXP(-XN+ 10 . 8/XKT I I ) ) + XNN ( I ) *6 T 51 

1 .4E-17*EXP(-XN+12 .O/XKT 1 1 ) ) ) / I.4.Q+1Q. 0*EXP I -2.38/XKTI I ) )+6. 0*EXP( - T 52 
23.57/XKTI I) ) ) +EQ5+EQ7) +BQ6 T 53 

GO TO 17 T 54 

13 IF (HV— 14.31 14,14,16 T 55 

14 X APNU I L » I ) S EPC*( EQ4+I 5. 16E— l 7*XNN ( I ) *EXPI -XN+10.8/XKT II) ) + XNN (11*6 T 56 

1* 4E— 17*EXPI— XN+12 • O/XKT 1 1 ) ) )/( 4.0+10. 0*EXPl-2.38/XKT (I) )+6.0*EXP<- T 57 
23.57/XKT( I) ) ) +EQ9+EQ7 ) T 58 

IF IHV-14.2) 15,15,17 T 59 

15 XAPNU(L,I)=XAPNU(L, II+BQ6 T 60 

GO TO 17 T 61 

16 XAPNUI L , I )=EPC*( EQ4+I 5. 16E— 17*XNN( I )*E XPI-XN+10.8/XKT (I ) )+XNNII)*6 T 62 

1.4E-17*EXPI-XN+12.0/XKTII) )+3.16E-17*XNNII> )/l4.0+10.0*EXP(-2.38/X T 63 
2KT1 II ) +6. 0*EXPI— 3.57/XKTI I ) ) )+EQ9+EQ7) T 64 

17 CONTINUE T 65 

RETURN T 66 

END T 67- 
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00 -Sj 0 s Ui 


SUBROUTINE FEMP (IHELP,IOPT) U 1 

C U 2 

COMMON /AAA/ A(6, 6 ) , A12 { 16 , 12 ) , A22 ( 7, 7 ) , ALP( 6 } , ALPT C 1 6 , 5) , BMT( 16) , U 3 
1C(6,16),CH(16,2),CP(16) ,DEGI(6) ,H( 1 6) , JAT( 1 6 , 5 I , JPH ( 1 6 ) ,KAT ( 6 ) , KQO U A 
2E ( 1 6 ) , R A< 16,2 ),RBO< 16,2), RCI 16,2), RDI 1 6 , 2 ) , RD1( 16, 2 > , REI 16, 2 ) ,RE1< U 5 
316, 2) , SOI 16J ,TB<3> ,TU( 16,2) ,TU2(16,2) , VN( 17 ) , VNE 1 15 ) , VNTC 17 ) , VNU< 1 U 6 
A6,6),W26(6),W3(16),Y(16),RF(16,2),RCi(16,2),HS, IG, IGMS, IGMSP, IGP, I U 7 
50N, IS, I SP ,1 SPNGP, I SPNG2 »N, NG» NP, PRESS , W27 , RHO, WM, SYU, ISPNG U 8 

COMMON /RAO/ YY ( 100 ) , TEE! 1 GO ) , FHV ( 20) , NHV, N Y,C2 , 1 Y U 9 

COMMON /RAD/ XNN(7,100) U 10 

COMMON /RAD/ NI HVC, FHVC ( 50 ) , AHV ( 50 ) , AH VL< 20 ) U 11 

COMMON /RAO/ C1,C3,C4,FLG,C5,FLG1 U 12 

COMMON /RAO/ YOELT , DELTD,FL1»FL2 U 13 

COMMON /RAD/ GEE! 8 ) , EP S (8 ) , NU( 20) , NDI 70 ) , FF ( 70) , GAMP { 70 » , WOL I 20 ) , F U 14 
1HVMI20) »FHVP(20) U 15 

COMMON /CIONCL/ F < 1 00 , 1 0) , F2 < 100, 10 ) , H VL( 70 ) , EP , K2, K1 , I FL , I YCON, I Q U 16 
11 »WMI »BIJI 100, 10) »GMIN( 100 ,10),GPLU(100,10),IAED U 17 

DIMENSION NICNI100) U 18 

DIMENSION PRESUOO) U 19 

DIMENSION HHI 100) U 20 

DIMENSION TGUSI15), DIGI15) U 22 

DIMENSION FFFI3) U 23 

I HELP=0 U 24 

W27=0.0 U 25 

DO 1 J=l, IS U 26 

1 W27=W27+ALP I J )*W26( J) U 27 

IF (IOPT.EQ.il GO TO 2 U 28 

HS=HS*W27 U 29 

2 CONTINUE U 30 

DO 8 1=1, N U 31 

K=1 U 32 

DO 5 KS= 1 » I S U 33 

IF (JATII ,K)-KAT(KS)) 6,3,4 U 34 

3 C(KS,I )=ALPT( I,K) U 35 

K=K+1 U 36 

GO TO 5 U 37 

C I KS, I ) =0 .0 U 38 

CONTINUE U 39 

GO TO 8 U 40 

DO 7 LS=KS, IS U 41 

C(LS,I)=0.0 U 42 

CONTINUE U 43 
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DO 13 1=1, N U 44 

DO 11 J=1 , 2 U 45 

IF (KODEm-1) 9,9,10 U 46 

9 CH( I, J)=RA( I, J>+RBO(IiJ) U 47 

GO TO 11 U 48 

10 CHI I, J )=RA( 1, J) U 49 

11 CONTINUE U 50 

W3( I )=0»0 U 51 

DO 12 L S= 1 , IS U 52 

12 W3( II=W3(I)+C(LS, I)*W26(LS) U 53 

13 CONTINUE U 54 

I G=N U 55 

IF (ION) 18,18,14 U 56 

14 DO 15 1=1, IG U 57 

15 VNE(I)=0.0 U 58 

DO 17 1=2, IG U 59 

DO 16 KS=2 » I S U 60 

16 VNE(I)=VNE(I)+C(KS»I)#DEGI(KS) U 61 

17 VNE ( I )=VNE I I)+C(l , I )*DEGI( 1 ) U 62 

18 NP=N+1 U 63 

ISP=IS+1 U 64 

IGMS=IG— IS U 65 

I GMSP= I GMS+1 U 66 

NG=0 U 67 

I SPNG= I S+NG U 68 

I SPNGP= I S PNG+1 U 69 

ISPNG2= ISPNGP+1 U 70 

IGP=IG+1 U 71 

DO 19 1=1, N U 72 

19 JPH( I )=0 U 73 

DO 20 1=1, IS U 74 

DO 20 K=l, IS U 75 

J=I G- I S+K U 76 

20 A(K, I) = C( I, J) U 77 

CALL INVERT (IS, A) U 78 

DO 22 K=l » IGMS U 79 

DO 22 J=1,IS U 80 

VNU(K»J)=0,0 U 81 

DO 21 1=1, IS U 82 

21 VNU ( K» J )=VNU( K, J ) +A ( I , J )*C ( I »K) U 83 

22 CONTINUE U 84 

DO 23 1=1, NP U 85 

23 VNT ( I )=0. U 86 
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VN(NP)=AMAXU VNT(NP)t0.02) 


u 

87 


00 24 1*1, N 


u 

88 


VN ( n = VNT(NP)*VNTU 1 


u 

89 

24 

VNC I |=AMAX1(VN{ I ) , 1.0E-06) 


u 

90 


00 27 1=1 , I GMS 


u 

91 


DO 25 J=1,IS 


u 

92 

25 

A12(I,J)=-VNU( I ? J ) 


u 

93 


DO 26 K=ISP,ISPNGP 


u 

94 

26 

A 12 ( I , K 1=0 ®0 


u 

95 

27 

CONTINUE 


u 

96 


DO 28 1=1 SP , I SPNGP 


u 

97 


DO 28 J= I SP , I SPNGP 


u 

98 

28 

A22 { I , J 1=0*0 


u 

99 


I A= 1 


u 

100 


J A= 1 


u 

101 


I X= 1 


u 

102 

29 

TBJA=TB( JAI 


u 

103 


CALL EQUIL (T BJA, PRESS , IX, N2P, .FALSE . ) 


u 

104 


TB { JA I =TB JA 


u 

105 


IF C IX— 1 1 30, 31,30 


u 

106 

30 

WRITE { 6, 60 ) TBJA 


u 

1C7 

31 

VA=G* 0 


u 

108 


DO 32 1*1, N 


u 

109 

32 

VA=VA+HU)*VNU > 


u 

110 


HSB = VA/ VN { NP J 


u 

111 


GO TO (33,34) , IOPT 


u 

112 

33 

HS=HSB/W27 


u 

113 


FN SN 36 


T 

A A 3 

34 

FFF ( JA1=HSB-HS 

i. 

u 

115 


T GUS ( I A ) = T8 JA 


u 

116 


D IG ( I A ) =F FF ( J A) 


u 

117 


IF (JA-2) 35,36,36 


u 

118 

35 

JA= JA+1 


u 

119 


GO TO 29 


u 

120 

36 

IF (ABS(FFF(2)/HS)~l.QE-04) 47,47,37 


u 

121 

37 

IF ( FFF ( 1 1— FFF ( 2 ) I 38,39,38 


u 

122 

38 

TP=TB(1)*(FFF(1)/ (FFF(1 l-FFF ( 2 ))) * C TB ( 2 >~TB ( 1 ) ) 


u 

123 


GO TO 40 


u 

124 

39 

TP=TB ( 1 ) 


u 

125 

40 

IF € TP 1 41,41,42 


u 

126 

41 

TP=0«75*TBI 1) 


u 

127 

42 

FFF S 1 ) = FF F { 2) 


u 

128 


TB( 1)=T8(2) 


u 

129 
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TB( 2)=TP U 130 

IF ( I A— 10 I 46,46,43 U 131 

43 WRITE (6,61) HS , C TGUS < LK) , D IG< LK) , LK= 1 , I A ) U 132 

I MIN=1 U 133 

DO 45 I X= 2, 11 0 134 

IF (ABS(DIG( IX) )-ABS(DIG( I MIN ) ) ) 44,44,45 U 135 

44 I MI N= I X U 136 

45 CONTINUE U 137 

TB(2)=TGUS( IMIN) U 138 

GO TO 47 U 139 

46 I A= IA +1 U 140 

GO TO 29 U 141 

47 TB2=TB{ 2) U 142 

CALL EQUIL ( TB2 ,PRE SS, I X»N2P , .TRUE . ) U 143 

TB( 2)=TB2 U 144 

VA=0« 0 U 145 

VB=0. 0 U 146 

VC=0»0 U 147 

V D= 0 . 0 U 148 

VE=0.0 U 149 

VF=1.0/TB(2)**2 U 150 

DO 57 1=1, N U 151 

J=1 U 152 

IF (TU( I , 1 )-TB( 2) ) 48,49,49 U 153 

48 J=2 U 154 

49 IF (KODE(I)-l) 50,53,50 U 155 

50 J=1 U 156 

IF (TU2 ( I , J )-TB(2 ) ) 51,52,52 U 157 

51 J-2 u 158 

52 CP( I)=1.98726*(R01( I , J ) *TB ( 2 ) -2 .0*REH I,J)*VF) U 159 

GO TO 54 U 160 

53 CP( I)=RC( I » J) +RD( I * J) *TB(2 )+RE( I, J )*VF U 161 

54 VA=VA+CP( I > *VN( I) U 162 

VC=VC+VN( I ) U 163 

IF (I-IG) 56,55,56 U 164 

55 V 8 =VC U 165 

56 VG= VN ( I )*W3( I ) U 166 

57 VD= VD+VG U 167 

WM= VD/ ( VB+VE/ ( 82. 0597*TB( 2 ) ) ) U 168 

RHO= (PRESS*WM)/(82.0597*TB(2) ) U 169 

VNT (NP )=VN( NP ) u 170 

DO 58 1=1, N U 171 

58 VNT (I ) = VN( I)/VN(NP) U 172 
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SYU=0.0 U 173 

00 59 M=1,N U 174 

BMT ( M ) = ( VNT (M ) *W3 (M))/W27 U 175 

59 SYU=SYU+BMT(MI U 176 

RETURN U 177 

U 178 
U 179 

60 FORMAT C21H0N-R DIO NOT CONVERGE, E15. 5) U 180 

61 FORMAT (19H0T DID NOT CONVERGE, E15. 5/ / ( 2E15 . 5 ) I U 181 

END U 182- 
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SUBROUTINE OEFIOJ ( N, I S ,C, VN, JPH, I OJ ) V 1 

DIMENSION C 16 ,16) » VNI17), JPHI16), I0JI13), X(20), CTI20), V(16), V 2 

1 TAU<7,7>, P(20,2O), ILU6), IPH(l6) V 3 

EQUIVALENCE <X,CT) V 4 

LN=N V 5 

LS= IS V 6 

LNP=LN+ 1 V 7 

LSP=LS+1 V 8 

DO 1 1=1, LN V 9 

1 VIIHN(t) V 10 

CALL ORDERV <LN,V,IPH) V 11 

M=LN/2 V 12 

DO 2 J= 1 » M V 13 

K=L N+l — J V 14 

II=IPHIJ> V 15 

IPH( J)=IPH<K) V 16 

2 IPHI K )= I I V 17 

L=l V 18 

DO 10 K=1 t LN V 19 

I ND= IPH I K ) V 20 

IF UPHIIND)) 3,3,10 V 21 

3 DO 4 1=1, LS V 22 

4 TAUIL, I >=C(I, IND) V 23 

DO 6 1=1, L V 24 

TP=0.0 V 25 

DO 5 J= 1 , LS V 26 

5 TP=TP + TAU(I ,J)*TAUIL, J) V 27 

P(I,L)=TP V 28 

6 P(L,I)=TP V 29 

DO 7 1=1,20 V 30 

C T ( I ) = 0 V 31 

7 CONTINUE V 32 

CALL DTLNEQ ( L , P, CT , X, DET ) V 33 

IF (ABS(DET)— 1.0E-05) 10,10,8 V 34 

8 IL(L)=IND V 35 

IF (L-LS) 9,11,11 V 36 

9 L=L+1 V 37 

10 CONTINUE V 38 

11 CALL ORDERI (LS,IL) V 39 

I OJ ( LS P ) = LNP V 40 

DO 12 J=1 , LS V 41 

12 10 J (J)=1L(J) V 42 

RETURN V 43 

END V 44- 
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SUBROUTINE OROERI <N,IPH) W I 

DIMENSION IPH( 16) W 2 

1 K=0 VI 3 

DO 3 J=2» N W 4 

IF { IPHI J-l J-IPHI J>) 3,3,2 W 5 

2 I T= I PH I J— 1 ) W fa 

I PH ( J-l 1 = IPH( J ) W 7 

I PHI J) = IT W 8 

K=1 W 9 

3 CONTINUE W 10 

IF IK) 4,4,1 W 11 

4 RETURN W 12 

END W 13- 
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SUBROUTINE ORDERV (N,V,IPH) X 1 

DIMENSION VU6), IPH( 16 ) X 2 

LN=N X 3 

DO 1 1=1, LN X 4 

1 IPH( I )= I X 5 

2 K=0 X 6 

DO 4 J=2» LN X 7 

IF (V( J-D-VC J)l 4,4,3 X 8 

3 VT=V( J-I ) X 9 

V(J-1»=V(J) X 10 

V(J»=VT X II 

I T= IPH ( J-l I X 12 

IPHC J— 11= IPH( J) X 13 

IPH( J)=IT X 14 

K=1 X 15 

4 CONTINUE X 16 

IF (K> 5,5,2 X 17 

5 RETURN X 18 

END X 19- 
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SUBROUTINE INVERT <N,A) 

DIMENSION AC 6 ,6 } , ZMC6,6I, N7C6J, 
EQUIVALENCE (N9Cl),N7m> 

DO 3 1= 1 , N 
DO 3 J=1 , N 
IF (I-J) 2,1,2 

1 ZM( I * J)=l .0 
GO TO 3 

2 ZM{ I , J >=0 *0 

3 CONTINUE 
N7( 1 )=0 

DO 15 1=1, N 

Z3=G.O 

DO 6 J=1 , N 

DO 4 K=1 , I 

IF (J-N7CK)) 4,6,4 

4 CONTINUE 
V=ABSCACJ,m 
IF C V- Z 3 1 6,6,5 

5 Z3= V 
N9C I)=J 
L=J 

6 CONTINUE 

IF CZ3-1.OE-07) 7,7,8 

7 Z3=l • 0E-07 

8 V=A(L,II 

DO 10 J=1 ,N 
ZMIL, J)=ZMCL, J)/V 
IF (J-I) 10,10,9 

9 A(L,J)=A(L,JI/V 

10 CONTINUE 

DO 14 J=1,N 
IF CJ-L) 11,14,11 

11 V=-A(J,I> 

DO 13 K=1,N 

ZM(J,K)=ZMCJ, K) +V*ZM( L , K) 

IF C K — I > 13,13,12 

12 A(J,K)=A(J,K) +V*A(L,K) 

13 CONTINUE 

14 CONTINUE 

15 CONTINUE 

DO 16 1=1 , N 

J=N9( I ) 

DO 16 K=1,N 

16 ACI ,K)=ZM( J,K) 

RETURN 

END 


Y 1 

N9C6) Y 2 

Y 3 

Y 4 

Y 5 

Y 6 

Y 7 

Y 8 

Y 9 

Y 10 

Y 11 

Y 12 

Y 13 

Y 14 

Y 15 

Y 16 

Y 17 

Y 18 

Y 19 

Y 20 

Y 21 

Y 22 

Y 23 

Y 24 

Y 25 

Y 26 

Y 27 

Y 28 

Y 29 

Y 30 

Y 31 

Y 32 

Y 33 

Y 34 

Y 35 

Y 36 

Y 37 

Y 38 

Y 39 

Y 40 

Y 41 

Y 42 

Y 43 

Y 44 

Y 45 

Y 46 

Y 47 
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1 

2 


SUBROUTINE EQUIL (TE,PR ,IX ,N2P, DOCB) Z 

LOGICAL DOCB Z 

Z 

Z 


COMMON /AAA/ A (6, 6) , AI2 ( 16 , 12 > * A22 C 7, 7 ) , ALP( 6) , ALPT i 1 6, 5) * BMT ( 16 1 1 Z 
1C <6 , 16 ), CHU6 , 2 ), CP (16), DEG K 6) ,H< 16) , JAT( 16, 5) , JPH ( 16) ,KAT(6),K0D Z 
2E<16) t RA( 16,2),RB0(16,2),RC(16,2),RD( 16 ,2 ) , RDH 16, 2 ) , RE( 16 , 2 ) ,RE1 ( Z 
316,2), SO( 16),T8<3 >,TU( 16, 2 ) , TU2( 16, 2) , VN( 17 ) , VNE< 16 ) , VNT< 1 7 ) ,VNU< 1 Z 
46,6>,W26(6>,W3( 16), YC 1 6 ) , RF (16,2) , RC1 ( 1 6,2 ) , HS , IG, I GMS, IGMSP, IGP, I Z 
50N,IS,ISP,I SPNGP, ISPNG2,N, NG, NP, PRESS, «27,RH0,WM,SYU, ISPNG Z 

COMMON /RAD/ YY(IOO) , TEE( 100) ,FHV( 20) , NHV,NY ,C2 , I Y Z 

COMMON /RAD/ XNN(7,100) Z 

COMMON /RAD/ NIHVC,FHVC (50 ) , AHVC 50 ) , AHVLI 20 ) Z 

COMMON /RAD/ Cl ,C3,C4* FLG, C5, FLG1 Z 

COMMON /RAD/ YDELT , DELTD,F L 1 ,FL2 Z 

COMMON /RAD/ GEE{ 8 > , EPS ( 8 ) , NU( 20) , ND< 70),FF(70) ,GAMP C 70) ,WOL( 20) , F Z 
1HVM (20),FHVP( 20) Z 

COMMON /CtONCL/ F ( 100 , 10) , F2C 100, 10 ) , HVL< 70 ) , EP , K2, K1 , l FL , I YCON, I Q Z 
il,WMI,8IJ( 100,10) , GMINC 100 , 10 ) ,GPLU{ 100 ,10 ) , I AED Z 

DIMENSION NICN(IOO) Z 

DIMENSION PRESUOO) Z 

DIMENSION HH( 100) Z 

DIMENSION DMUC16), VLNKU6), STEP(17), EZ(17), AA<6), A21(7,17), A Z 

17(20,20), B( 20) , I0JI13), XB ( 20 ) , XX( 20 ) , DY{17) Z 

EQUIVALENCE (XB,XX) Z 

KCBAL=0 Z 

IF (DOCB) GO TO 54 Z 

V ANALOG ( TE/3000 * ) Z 

VB=TE-3000. Z 

VC- ( TE+3000 • ) /2 #0 Z 

V D= 1 * 0 / ( T E *30 00 * ) Z 

VE=VC*VD**2 Z 

JC=1 Z 

JD= 1 Z 

VF-G *0 Z 

KERRS V=KERR Z 

KERR=1 Z 

DO 1 1=1, IG Z 

VF=VF + VN( I ) Z 

IF (VF-PR) 2,4,2 Z 

VF=PR/VF Z 

DO 3 1=1, NP Z 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 
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3 VN( I )=VF*VN( I ) Z 44 

4 DO 5 1 = 1, IG Z 45 

5 Y ( I ) = A L OG ( VN I I)) Z 46 

DO 14 1=1, N Z 47 

J=1 Z 48 

IF (TE-TUI I»J ) ) 7,7,6 Z 49 

6 J=2 Z 50 

7 IF (KODEI 1 1—1 ) 8,8,9 Z 51 

8 HIIKHII ,JJ+V8*(RC(l, J)+R0( I, J>*VC+RE( I,J)*VD> Z 52 

SD( I) = RF( I,J>+RCC I, J)*VA+VB*(RD(I, J H-RE 1 1 , J ) *VE ) Z 53 

DMU( I)=H( I )-TE*SD( I) Z 54 

GO TO 14 Z 55 

9 J=1 Z 56 

IF I TE-TUI I » J ) I 11,11,10 Z 57 

10 J=2 Z 58 

11 DMU(I)=CH(I,J)-1.98726*{RC(I,J)*TE+RD( I ,J»*TE**2+RE( I, J)/TE) Z 59 

J=1 Z 60 

IF ( TE-TU2 ( I , J) > 13,13,12 Z 61 

12 J=2 Z 62 

13 SD( II=1.98726*(RC1( I , J I+RD1 ( I , J )*TE+RE 1 C I , J ) /TE**2) Z 63 

HI I)=DMU( n+TE*SD(I) Z 64 

14 CONTINUE Z 65 

VA=-l,0/t 1.98726+TE) Z 66 

DO 16 1=1 , IGMS Z 67 

VB=0.0 Z 68 

00 15 LS=1» IS Z 69 

J=IGMS+LS Z 70 

15 VB=VB+VNUl I,LS)*DMU(J» Z 71 

16 VLNKI n=VA*IOMU(I)-VB) Z 72 

CALL OEFIOJ (N,IS,C,VN, JPH,IOJ) Z 73 

CALL MASS AL ( VN,I OJ,C, ALP, JC,0 .9, IS,N,PR» IG) Z 74 

17 EL=0.0 Z 75 

ENL=0.0 Z 76 

DO 18 I=1,NP Z 77 

18 STEP(I)=1.0 Z 78 

00 19 1=1, IG Z 79 

19 Y(n=ALOGlVNU») Z 80 

DO 23 1=1, IGMS Z 81 

VA=0.0 Z 82 

DO 20 L S= 1 » I S Z 83 

J=I GMS+LS Z 84 

20 VA=VA+VNU( I,LSJ*Y( JI Z 85 

EZ( I)=VLNK< I)-YU)+VA Z 86 
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A12(I, ISPNG2)=EZ( I> Z 87 

IF (VN( D-1.0E-30) 21,21,22 Z 88 

21 STEP(1 ) =0 .0 Z 89 

EL* AMAX1( EL,EZ(I)) Z 90 

GO TO 23 Z 91 

22 EL=AMAX1(EL»ABS(EZ< ill) Z 92 

23 CONTINUE Z 93 

DO 24 1=1, IS Z 94 

AA(I)=0.0 Z 95 

DO 24 J=1 * IGMS Z 96 

A21(I,J)=C( I, J>*VN< J> Z 97 

24 AA( n = AA(I)+A2HI,J) Z 98 

DO 25 J=1 , IGMS Z 99 

25 A2K ISPNGP, Jl=VN( Jl Z 100 

DO 26 1=1, IS Z 101 

DO 26 J=1 , ISPNG Z 102 

K=I GMS+ J Z 103 

A22 (I,J l=C( I,KJ*VN(K) Z 104 

26 AAl I )= AA( I ) +A22CI , J I Z 105 

DO 27 J=1,IS Z 106 

K=IGMS+J Z 107 

27 A22 ( ISPNGP, J)*VN( K) Z 108 

00 28 1=1, IS Z 109 

A22 ( I , ISPNGP l=— AA ( I ) Z 110 

VA=( VN( NP)*ALP( I) )/AA( I ) Z 111 

K=IGMS+ 1 Z 112 

EZ(K1=AL0G(VA> Z 113 

A22 ( 1 , ISPNG2)=AA( I ) *EZ (K) Z 114 

ENL=AMAX1 (ENL , ABS t EZ{ K) ) ) Z 115 

28 CONTINUE Z 116 

PB=0.0 Z 117 

DO 29 1=1, IG Z 118 

29 PB=P8+VN( I ) Z 119 

EZ ( NP ) = ALOG(PR/PB J Z 120 

A22(ISPNGP, ISPNG2 )=P8*EZ(NPJ Z 121 

ENL* AM A XI { ENL , ABS (EZtNP ) ) ) Z 122 

IF (AMAXltEL, ENL )— 1 .0E-05 ) 54,54,30 Z 123 

30 00 32 1=1, ISPNGP Z 124 

DO 32 K=l, ISPNGP Z 125 

VA=0.0 Z 126 

DO 31 J=1 , IGMS Z 127 

31 VA=VA+A21t I,J»*A12(J,K>*STEPCJI Z 128 

32 AT( I,K)=A22( I,K)-VA Z 129 
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00 34 1=1 » I SPNGP 

Z 

130 


VA=0.0 

Z 

131 


00 33 J = 1 , I GM S 

Z 

132 

33 

VA= VA+ A21 4 I , J I *A1 2 4 J , I S PNG2 I 

z 

133 

34 

BCn = A22( If IS PNG2 )-VA 

z 

134 


CALL OTLNEQ 4 1 SPNGP, AT , 8, XB, VA) 

z 

135 


XX4 ISPNG2 )=-l . 

z 

136 


VB=G*0 

z 

137 


DO 37 1=1 t I SPNGP 

z 

138 


J= I GMS + I 

z 

139 


DY( J)=XX4 n 

z 

140 


IF (I-IS) 36f 36t3 5 

z 

141 

35 

IF (STEP4JI) 37,37, 36 

z 

142 

36 

VB=AMAX14 VB,ABS(XX( I) I I 

z 

143 

37 

CONTINUE 

z 

144 


DO 41 1 = 1 * IGMS 

z 

145 


VA=0.0 

z 

146 


DO 38 J = 1 * I SPNG2 

z 

147 

38 

VA=VA+A124I, JI*XX4JI 

z 

148 


DY ( I )=- VA 

z 

149 


IF C ST EP( I ) ) 39» 39f 40 

z 

150 

39 

DY( 1 1 = A 12 ( I f I SPNG 2) 

z 

151 


GO TO 41 

z 

152 

40 

V 8= AM AX! 4 A8S4 VA)t VB) 

z 

153 

41 

CONTINUE 

z 

154 


IF 4 VB-2* 3052 1 44,44,42 

z 

155 

42 

VA=2.3052/VB 

z 

156 


DO 43 1=1 tNP 

z 

157 

43 

DY4 I) = VA*DY(n 

z 

158 

44 

DO 45 1=1, IG 

z 

159 


VN4 1) = AMAX14VN( I) *£XP 4 OY 41 ) ) , 1 #Q£-30) 

z 

160 

45 

CONTINUE 

z 

161 


DO 46 I=IGP,NP 

z 

162 

46 

VN4I)=AMAX1(VN( I)*EXP(DY(I )),5.GE-10) 

z 

163 


IF ( JC— 50 ) 48,48,47 

z 

164 

47 

I X=2 

z 

165 


RETURN 

z 

166 

48 

IF 4 J D- 3 ) 50,49,49 

z 

167 

49 

CALL DEFIOJ 4 N, IS ,C , VN, JPH, IOJ) 

z 

168 


J D= 1 

z 

169 


GO TO 51 

z 

170 

50 

JD= JD+1 

z 

171 

51 

IF 4(1*111* END— ELI 53,52,52 

z 

172 
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52 CALL MASBAL { VN , I OJ,C, ALP, JC,0*9, I S, N , PR, I G) Z 173 

53 JC=JC+1 Z 174 

GO TO 17 Z 175 

54 IF UNOT.DGCB) RETURN Z 176 

IF ( ION.EQ.Q.OR,IX*EQ.2) RETURN Z 177 

VNSAVE^O.O Z 178 

00 55 M=2,N Z 179 

55 VNSAVE=VNSAVE+VN<M)*VN£<MI Z 180 

CHGERR=VNSAVE/VN(1)-1.0 Z 181 

IF <ABS(CHGERRH1.0E-04) 56,56,57 Z 182 

56 RETURN Z 183 

57 KCB AL-KCBAL+l Z 184 

IF (KCB AL *GT. 3 ) RETURN Z 185 

CALL CHGBAL (TE, PR, KRAP,OMU) Z 186 

IF (KRAP* EQ*2 ) RETURN Z 187 

JC= 1 Z 188 

GO TO 17 Z 189 

END Z 190- 
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1 


2 

3 


4 


5 


6 


SUBROUTINE CHG8AL ( TE, PR,KRAP»DMU) AA 1 

REAL KSU8 P,LGGKP AA 2 

DIMENSION ALPTC (16,7), JATC(16,7), DMUCIC16), VNECB<16), VNMN(16), AA 3 

1 DMUCN(16), KSUBP ( 16) , AK(1G), INDEX<16), VNOLD ( 16) , DMUC 16 ) AA 4 

COMMON /AAA/ A(6 f 6) , A12 ( 16, 12 ) , A22 ( 7, 7) , ALPC 6) , ALPTi 1 6 , 5) , BMT{ 16) * AA 5 
1C (6,161, CH( 16, 2), CP (16) »DEGI<6) » H ( 161 , J ATU6 , 5 J , JPH ( 16 ) * KAT< 6) ,KGD AA 6 
2E(16),RAC 16,21, R8 0(16, 2), RC(16,2),RD( 16,2 ) , RDH 16,2) , RE ( 16 , 2 1 ,RE1 ( AA 7 
316,2),SD( 16),TB(3),TU<16,2) ,TU2(16,2i ,VN(17) , VNEU6 ) , VNTU7 ) , VNU ( 1 AA 8 
46f6)fW26(6),W3(16),y(16),RF(16,2),RCl (16,2) , HS, IG, I GMS, IGMSP , IGP , I A A 9 


5QN, IS , I SP , I SPNGP, IS PNG 2 ,N, NG,NP, PRESS ,W27, RHO, WM, SYU, ISPNG A A 10 

COMMON /RAD/ YYClGOi, T££( 100 >, FHVI20 ), NHV, NY, C2, I Y AA 11 

COMMON /RAD/ XNN(7,100) AA 12 

COMMON /RAO/ NIHVC, FHVC < 50 ) , AHV { 50 i , AHVL ( 20 ) AA 13 

COMMON /RAD/ Cl ,C3, C4,FLG, C 5,FLG1 AA 14 

COMMON /RAD/ YDELT, DELTD,F(1 1,FL2 AA 15 

COMMON /RAD/ GEE( 8) ,EPS(8 > , NU(20) ,ND( 70) , FF < 70) ,GAMP( 70 ) , tfGL< 20) , F AA 16 

1HVM(20) ,FHVP(20) AA 17 

COMMON /CIONCL/ F ( 100, 10) , F2( 100, 10 ) , HVL( 70 ) , EP , K2, K1 , I FL , I YCON, I Q AA 18 

11* WMI,BIJ( 100,10) , GMIN( 100 , 10) ,GPLU( 100, 10) , I AED AA 19 

DIMENSION NICN(IOO) AA 20 

DIMENSION PRES(IOO) AA 21 

DIMENSION HH(IOO) AA 22 

PB=G* 0 AA 24 

DO 1 1*1, IG AA 25 

PB=PB+VN(I) AA 26 

J=1 AA 27 

DO 4 1=2, IG AA 28 

IF (ABS(VNE( I ) )— 1 *GE-20 ) 4,4,2 AA 29 

DO 3 K=1 , 5 AA 30 

ALPTC( J,K)=ALPTU ,K) AA 31 

J ATC( J* K) = JAT ( I ,K ) AA 32 

DMUC I f J )=DMU( I) AA 33 

VNECBt J ) = VNE( El AA 34 

I NDEX f J 1= I AA 35 

J=J + 1 AA 36 

CONTINUE AA 37 

NQI0N=J-1 A A 38 

DO 5 J*l, NOION AA 39 

DMUCN( J)=Q*0 AA 40 

DO 9 1=2, IG AA 41 

IF (ABS(VNEd) )-l„0E-20) 6,6,9 AA 42 

DO 8 J=l, NOION AA 43 
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JIN=1 


AA 

44 


JNU=1 


AA 

45 


IF { JATU,1)*EQ*0 ) JNU=2 


AA 

46 


IF C JATCt JylKEQ.Ot JIN-2 


AA 

47 

7 

DIFF=ABS( ALFT( X, JNUJ-ALPTCi J, 

JIN) ) 

AA 

48 


IF ( < J AT C I * JNUI .NE.JATC(J,JIN)) .OR. (OIFF.GT • 1.0E-2D ) ) GO TO 8 

AA 

49 


JNU=JNU+1 


AA 

50 


J I N= J I N+l 


AA 

51 


IF (MAXOC JNU, JIN) *LT,6> GO TO 

7 

AA 

52 


VNMN( J) =VN( I I 


AA 

53 


DMUCNC JI=DMU( l) 


AA 

54 

8 

CONTINUE 


AA 

55 

9 

CONTINUE 


AA 

56 


MI S ING=0 


AA 

57 


00 10 J-l *NGI ON 


AA 

58 


IF (ABSCOMUCNI J)) .GT.UGE-10) 

GO TO 10 

AA 

59 


WISING- MISING+1 


AA 

60 


1=1 NDEX { J J 


AA 

61 

10 

CONTINUE 


AA 

62 


IF <MISING.EQ*0) GO TO 11 


AA 

63 


CALL DUMP 


AA 

64 

11 

RTR=— 1 *0/ ( TE* 1*98 726) 


AA 

65 


DO 12 J-l , NOI ON 


AA 

66 


LOGKP=RTR* (DMUCl ( J )-DMUCNC J )+VNECB( J) *DMU( 1 > ) 

AA 

67 


KSUBP( J)=EXP( LOGKP) 


AA 

68 

12 

1 = 1 NOEX ( J 1 


AA 

69 


DO 15 KK=1,5 


AA 

70 


AK(KK)=0.0 


AA 

71 


DO 14 J=l, NOION 


AA 

72 


KVNECB= IF IX ( VNECB I J )+0 .001 ) 


AA 

73 


IF (VNECB( JI.LT.O.O) KVNECB=IFI XI VNECB ( J ) -0 .001) 

AA 

74 


IF i ( KK-3 )+KVNECB ) 14,13,14 


AA 

75 

13 

AM KK)=AK(KK)+VNECB(J)*KSUBP( J)*VNMN( J ) 

AA 

76 

14 

CONTINUE 


AA 

77 

15 

CONTINUE 


AA 

78 


AK ( 4 ) = AK ( 4 ) — 1 * 0 


AA 

79 


XNEW-ABSCVNIl 1 ) 


AA 

80 


KRU0D=1 


AA 

81 

16 

DO 18 LL=1*50 


AA 

82 


L = LL 


AA 

83 


FOF X=0*Q 


AA 

84 


FPOF X- 0 ®0 


AA 

85 


X=XNEW 


AA 

86 


89 




DO 17 KK=1,5 

AA 

87 


VKK=KK-3 

AA 

88 


TEST1=VKK*AL0G<X) 

AA 

89 


TEST2=< VKK-1.0)*ALGG(X) 

AA 

90 


T ERM1- 1 #0 E+36 

AA 

91 


TERM2=l,0E+36 

AA 

92 


IF (TEST1*LT.83*Q) TERM1=X**VKK 

AA 

93 


IF (TEST2*LT*83*0) TERM2=X** < VKK-1 • 0) 

AA 

94 


FOFX=FOFX+AK(KK)*TERMl 

AA 

95 

17 

FPOFX=FPOFX+VKK*AK( KKI *TERM2 

AA 

96 


XNE W-X— ( FOFX/ FPOFX ) 

A A 

97 


IF (XNEW. LT. 0.0) XNEW^X/2. 0 

AA 

98 


ERRGR=1.0-XNEW/X 

AA 

99 


IF ( ABS (ERROR )— 1. 0E— 06 ) 20,20,18 

AA 

100 

18 

CONTINUE 

AA 

101 

19 

KRAP=2 

AA 

102 


RETURN 

AA 

103 

20 

IF { XNEW ) 21,22,22 

AA 

104 

21 

XNEW-AB S( XNEW) 

AA 

105 


KRUDD=KRUDD+1 

AA 

106 


IF (KRUDD.GT • 2 ) GO TO 19 

AA 

107 


GO TO 16 

AA 

108 

22 

X^XNEW 

AA 

109 


VN< 1 )=X 

AA 

110 


DO 23 J»1,N0IGN 

AA 

111 


1=1 NDEX { J ) 

AA 

112 


VNOLEM I )= VN( I ) 

AA 

113 


VN< I)=KSU8P(J)*VNMN( J) *X** (-VNECB { J I) 

AA 

114 

23 

VN( n=AMAXl(VN( I) ,9.999E-3G ) 

AA 

115 


KRAP^l 

AA 

116 


RETURN 

AA 

117 


END 

AA 

118- 
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sO 00 -g & \J\ > OJ N 


SUBROUTINE ZHV ( HV, ZO, ZN, Z I , ZC ) AB L 

X=HV AB 2 

IF IHV-9.82) 1,1,2 AB 3 

I ZO=C. 99997956-0. 31554804*X+2.8245479E-02*X**2+6.6773283E-03*X**3-3 AB 4 

1.6445854E-03*X**4+8.0580698E-04*X**5-7. 7086374E-05* X**6+2. 668133E- AB 5 
206*X**7 AB 6 

GO TO 3 AB 7 

Z0=(X/9.82>**3 AB 8 

IF (HV— 8.351 4,4,5 AB 9 

ZN= 1.00 0148—0. 418 35346*X+0. 16803591*X**2-9. 7794579E-02*X**3+3. 3546 AB 10 
1351E-02*X**4-5.6093 534E-03*X**5-«-4.515535E-04*X**6-1.4035845E-05*X* AB 11 
2*7 AB 12 

GO TO 6 AB 13 

ZN= ( X/8 .35 ) **3 AB 14 

X=HV/4 .0 AB 15 

IF (X-6.6J 7,7,8 AB 16 

Z 1=1. 0003794-0. 29 64 7668*X+7.5052416E-02*X**2-1. 702948 lE-02*X**3+3 . AB 17 
12795539E-03*X**4-2. 1284692 E—04*X** 5 AB 18 

GO TO 9 AB 19 

ZI= ( X/ 7 .37 ) **3 AB 20 

X=H V AB 21 

IF (X-7.37) 10,10,11 AB 22 

10 ZC=0. 99743674-0. 43418122*X+8.5313141E-02*X**2-1.3939168E-02*X**3+4 AB 23 

1 ,0385449E-03*X**4— 5.4264246E-04*X**5+2. 812126 IE— 05*X**6-3. 8 83 5 29 8 E AB 24 
2-07*X**7 AB 25 

GO TO 12 AB 26 

II ZC=(X/7.37)**3 AB 27 

12 RETURN AB 28 

END AB 29- 
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SUBROUTINE LTP UTSW) 


AC 

1 



AC 

2 

COMMON /AAA/ A(6,6) , A12 U6 , 12) , A22 < 7 , 7) , ALP ( 6 > , ALPT ( 1 6 , 

5) , BMT ( 16 ) , 

AC 

3 

lC(6»16)»CH(16»2)»CP(16)tDEGI(6)»H( 16) , JAT<16,5) ,JPH(16) 

,KAT ( 6 ) »KQO 

AC 

4 

2E 1 16), RAC 16, 2 ),RBO( 16,2 ),RC( 16,2) ,RQ( 16,2) ,RD1( 16,2) ,RE(16,2),RE1( 

AC 

5 

316,2),SD(16)»TB(3).TU( 16, 2 ) , TU2 ( 16, 2 ) ,VN< 17 ) , VNEC 16 ) , VNT( 1 7 ) , VNUC 1 

AC 

6 

46,6 ) ,W26(6) ,W3( 16) ,Y(16) ,RF(16,2) ,RC1 ( 16,2) ,HS, IG, IGMS, 

IGMSP, IGP, I 

AC 

7 

50N, IS, ISP, I SPNGP, ISPNG2, N»NG,NP, PRESS, W27,RH0,WM,SYU, ISPNG 

AC 

8 

COMMON /RAD/ YY(100),TEE(1QQ),FHV(2O>, NHV, NY , C2, I Y 


AC 

9 

COMMON /RAD/ XNN(7,100> 


AC 

10 

COMMON /RAC/ NIHVC, FHVC (50 ) * AHV (50 ) * AHVL( 20 ) 


AC 

11 

COMMON /RAD/ Cl ,C 3, C4, F LG, C 5 , FLG1 


AC 

12 

COMMON /RAD/ YOELT , DELTD, F L 1, FL2 


AC 

13 

COMMON /RAO/ GEEf 8 ) , EPS ( 8 1 , NU( 20) , ND{ 70 ) , FF { 70) * GAMP ( 70 

) » WOL ( 20 ) , F 

AC 

14 

1HVM C 20 ) * FHVP ( 20 ) 


AC 

13 

COMMON /CIONCL/ F < 100* 10) , F2 (100, 10 ), HVLC 70 ) f EP »K2, K1 , I 

FL, IYCQN, IQ 

AC 

16 

U*WMI,BIJdOO f 10> fGMINIl 00,10) , GPLLM 100,10),! AED 


AC 

17 

DIMENSION NIC N (100) 


AC 

18 

DIMENSION PRESdOOl 


AC 

19 

DIMENSION HHdOO) 


AC 

20 

ALP (1 ) =5* 5 


AC 

22 

ALP ( 2 ) = 1 ® 4375 


AC 

23 

I T S &= 1 


AC 

24 

ION=0 


AC 

25 

I S=2 


AC 

26 

N=4 


AC 

27 

W26 ( 11= 14*008 


AC 

28 

W26 { 2 ) = 16 *0 


AC 

29 

KAT ( 11=7 


AC 

30 

K AT C 2 1 — 8 


AC 

31 

D£GI(n=0,0 


AC 

32 

DEG 1(2 )=0 «Q 


AC 

33 

KODEC 1 ) -1 


AC 

34 

DO 1 J=2,4 


AC 

35 

KODEI J 1 =2 


AC 

36 

DO 2 J=l,4 


AC 

37 

DO 2 K=l,5 


AC 

38 

ALPTI J , K) =0 *0 


AC 

39 

JATC J,KS=0 


AC 

40 

ALPT (1,11—2 »0 


AC 

41 

JAT (1,1 1=8 


AC 

42 

ALPTC2, 11=1*0 


AC 

43 
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JAT(2»11=j 

AC 

44 

ALPT(3, 11=1.0 

AC 

45 

JAT (3*1 1=8 

AC 

46 

ALPT(4, 11=2.0 

AC 

47 

JATC4,1)=7 

AC 

48 

RA( 1,11=0.0 

AC 

49 

R BO (1,11=2.551 97E+04 

AC 

50 

RC(1, 11=8. 12927 

AC 

51 

RD(1,11=4. 85368E-04 

AC 

52 

RE( 1,1 )=-2. 75743E+05 

AC 

53 

RF(1, 11=67. 9798 

AC 

54 

TU( 1, 11=3000. 

AC 

55 

RA( 1,21=0.0 

AC 

56 

RB0(l,21=2.55239E+04 

AC 

57 

RC( 1,21=10.2348 

AC 

58 

RD( 1,21=4. 5738E-05 

AC 

59 

RE(1,21=— 7.73465E+06 

AC 

60 

RFI 1,21 =67.9813 

AC 

61 

TU(1, 21=6000. 

AC 

62 

RA(2,ll=0.1125E+06 

AC 

63 

RC(2, 11=19. 8515 

AC 

64 

RD( 2,1 1=0.650822E— 03 

AC 

65 

RE( 2, 1 l=-0. 156336E+07 

AC 

66 

TU( 2,11=5000. 

AC 

67 

R A( 2, 2 1 =0. 1 12 5E+06 

AC 

68 

RC(2, 21=23. 2779 

AC 

69 

R0(2,21=0. 190986E— 03 

AC 

70 

RE(2»2)=— 0.3310C3E+08 

AC 

71 

TU(2, 21=24000. 

AC 

72 

RC1(2, 11=22.4633 

AC 

73 

RD1 (2, 1 )=0.622485E— 03 

AC 

74 

REl(2,ll=-0. 165345E+07 

AC 

75 

TU2(2, 11=6000. 

AC 

76 

RC1(2, 21=26. 0356 

AC 

77 

RD1(2,21=0. 229541E-03 

AC 

78 

RE1 ( 2, 2 1=— 0.501 696E +08 

AC 

79 

TU2(2, 21=24000. 

AC 

80 

RA( 3, 1 1=0 . 589 85E+05 

AC 

81 

RC(3, 11=20. 8169 

AC 

82 

RD( 3, 1 1=0 .655447E— 03 

AC 

83 

RE ( 3,1 l=-0.162744E+07 

AC 

84 

TU( 3,11=5000. 

AC 

85 

RA( 3,21=0.58985E+05 

AC 

86 
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RC(3, 21=24.5103 

AC 

87 

RD(3»2)=0. 158277E— 03 

AC 

88 

RE(3»2)=— 0.351471E+08 

AC 

89 

TU( 3,2 >=24000. 

AC 

90 

RC1(3, 11=23.3381 

AC 

91 

R01 (3,1 )=0.659562E— 03 

AC 

92 

RE1 ( 3, 1 )=— 0.1 55587E+07 

AC 

93 

TU2(3, 11=5000. 

AC 

94 

RC1(3, 21=26. 8167 

AC 

95 

RD1 ( 3, 2 )=Q. 186999E-03 

AC 

96 

R El (3,2 )=-0.322302E+08 

AC 

97 

TU2(3, 21=24000. 

AC 

98 

RA(4, 11=0.0 

AC 

99 

RC ( 4, 1 1 =24. 91 73 

AC 

100 

RD(4,1)=1.08123E— 03 

AC 

101 

RE ( 4, 1 1 =-2.21 576E+06 

AC 

102 

TU( 4, 11=5000. 

AC 

103 

RA(4, 21=0.0 

AC 

104 

RC(4, 21=30. 6229 

AC 

105 

R0(4,21=3 .01129E— 04 

AC 

106 

RE ( 4,2 1 =- 5. 19808E+07 

AC 

107 

TU(4, 21=24000. 

AC 

108 

RC1(4, 11=28.2374 

AC 

109 

R01 (4, 1 )=1.32899E— 03 

AC 

110 

RE1 (4, ll=-2.1492E+06 

AC 

111 

TU2 (4, 11=4000. 

AC 

112 

RC1(4,2)=33.5174 

AC 

113 

RD1 (4, 2 1=4.39119E— 04 

AC 

114 

RE1 (4, 21=— 3.14737E+07 

AC 

115 

TU2(4, 21=24000. 

AC 

116 

RETURN 

AC 

117 

END 

AC 

118— 
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eg m 4* 


SUBROUTINE HAFACE (IFL,01,D2) 


AD 

1 

DIMENSION ALOWER! 40 ) , TABLE! 1) , 

ALINIT! 40) , X1STEP! 40 ) , XNU(40), 

X AD 

2 

1NUBAR (40) t SAVE(40), XSTAR(40>, 

DX ( 40 ) , DELT ( 40 ) 

AD 

3 

IF ! IFL.EQ.2) GO TO I 


AD 

4 

ARG5=D1 


AD 

5 

TF=D2 


AD 

6 

N0RD=1 


AD 

7 

NRT AB=0 


AD 

8 

DTP=(TF—ARG5) *.l 


AD 

9 

HMIN=DTP*.02 


AD 

10 

ARG4=DTP*.l 


AD 

11 

HMA X=DTP 


AD 

12 

RMAX=l.E-5 


AD 

13 

TP= ARG5 


AD 

14 

ALI NIT! 1 ) =0. 


AD 

15 

RETURN 


AD 

16 

CALL INOCAL ! XNU, ALINIT ,-l ARG5 

,NTAB) 

AD 

17 

NTAB=1 


AD 

18 

SEN S=0« 


AD 

19 

GO TO 36 


AD 

20 

SENSR=0. 


AD 

21 

IF (TF-ARG5) 44,4,4 


AD 

22 

DO 5 1= l » NORD 


AD 

23 

ALOWER! I ) =A'LI NIT! 11) 


AD 

24 

SAVE! I )=ALINIT( I ) 


AD 

25 

XNUBAR! I)=ALINIT(I) 


AD 

26 

ARG7=ARG4/2.G 


AD 

27 

ARG8=ARG7/2 .0 


AD 

28 

ARG6=ARG4 


AD 

29 

SW3=— 1 . 


AD 

30 

ARG9=ARG5 


AD 

31 

ARG10=ARG6 


AD 

32 

ARG11=ARG7 


AD 

33 

CALL INOCAL ( X STAR » XNU BAR » 0 .0» ARG9 » NTAB I 

AD 

34 

DO 7 I* I, NORD 


AD 

35 

DX( I) = XSTAR(II*ARG10 


AD 

36 

XNUBAR ( n = XNUBAR! I) + !DX(I)/2.0) 


AD 

37 

ARG9=ARG9+ARGli 


AD 

38 

CALL INOCAL ( XSTAR, XNUBAR, 0 .0, ARG9 , NT AB ) 

AD 

39 

DO 8 1=1, NORD 


AD 

40 

TEMP=XSTAR( I) +ARG10 


AD 

41 

DX( I )=DX{ I )+( TEMP+2.0) 


AD 

42 
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8 

XNUBARI U = SAVE( I) + C TEMP/2.0) 

AD 

43 


CALL INOCAL ( XSTAR, XNUBAR,0 .0, ARG9 , NTA8 ) 

AD 

44 


00 9 I -It NORD 

AD 

45 


TEMP=X$TAR { I ) *ARG10 

AD 

46 


0X{ I ) = DX ( E)+(TEMP*2.0) 

AD 

47 

9 

XNU8AR ( I ) =TEMP+SAVE ( I ) 

AD 

48 


ARG9=ARG9+ARG11 

AD 

49 


CALL INOCAL ( XSTARt XNUBAR ,0.0, ARG9 , NTA8 ) 

AD 

50 


DO 10 1 = 1 1 NORD 

AD 

51 

10 

XSTAR C I ) = C C tXSTARU)*ARG10)+DX( m/6.G)+SAVE{I) 

AD 

52 


IF ( SW3 ) 1 1 ? 13 1 15 

AD 

53 

11 

SW3=0,0 

AD 

54 


ARG9=ARG5 

AD 

55 


ARG10= ARG7 

AD 

56 


ARG 11=ARG8 

AD 

57 


00 12 1=1, NORD 

AD 

58 


XNU8AR(I)=SAVE{I) 

AD 

59 

12 

X1STEPC I)=XSTARU) 

AD 

60 


GO TO 6 

AD 

61 

13 

SW3=1 « 0 

AD 

62 


ARG9=ARG5+ARG7 

AD 

63 


DO 14 1=1, NORD 

AD 

64 


XNU8AR CI) = XSTARCI ) 

AD 

65 

,14 

SAVE 1 1 ) = XNUBARC 1 1 

AD 

66 


GO TO 6 

AD 

67 

15 

DO 16 1=1, NORD 

AD 

68 


DELTCI ) = < XSTAR( I ) — X1STEP( I ) 1/15.0 

AD 

69 


XNUC I l = XST AR I 1 )+DELT( I ) 

AD 

70 

16 

A LOWER ( I )=XNU( I ) . 

AD 

71 


QMAX=0 • 

AD 

72 


DO 25 1=1, NORD 

AO 

73 


IF CABSCXNUCI M-l.O) 17,18,18 

AD 

74 

17 

QUG=AB$ C DELTC U) 

AD 

75 


GO TO 19 

AD 

76 

18 

QUO= A8S ( DELTCI) /XNU C I)) 

AD 

77 

19 

IF i QM AX— QUO ) 20,20,21 

AD 

78 

20 

QMAX=QUO 

AD 

79 

21 

IF { RMAX— QUO) 22,22,25 

AD 

80 

22 

IF (SENS) 24,23,24 

AD 

81 

23 

ARG4=ARG4/2.0 

AD 

82 


IF CARG4-HMIN) 43,43,2 

AD 

83 

24 

QMAX=RMAX 

AD 

84 

25 

CONTINUE 

AD 

85 
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DO 26 1=1 t NGRD 

AD 

86 

26 

ALI NIT { n=ALOWER(Il 

AD 

87 


ARG5=ARG5+ARG4 

AD 

88 


IF (TP-ARG5) 31,31,27 

AD 

89 

27 

IF { SENSR ) 28,40,28 

AD 

90 

28 

IF (TABLE ( NTAB ) +DTP— TP ) 30,30,29 

AD 

91 

29 

TP=TABLE( NTAB J+DTP 

AD 

92 

30 

GO TO 32 

AD 

93 

31 

TP=TP+DTP 

AD 

94 

32 

CALL INOCAL ( XSTAR, XNU, 0,0 , ARG9 , NTAB) 

AD 

95 


SUB=1,0 

AD 

96 


CALL INOCAL ( XSTAR, XNU , SUB , ARG9 , NTAB) 

AD 

97 


IF (SUB-1 5,0) 40,44,40 

AD 

98 

33 

ARG4=HMAX 

AD 

99 

34 

IF (SENSR) 35,36,35 

AD 

100 

35 

NTAB=NTAB+1 

AD 

101 

36 

IF (NTAB-NRTAB) 37,37,2 

AD 

102 

37 

IF (TABLE (NTAB) -ARG5-ARG4) 38,39,2 

AD 

103 

38 

ARG4=TABLE(NTAB)-ARG5 

AD 

104 

39 

SENSR= 1 , 

AD 

105 


GO TO 3 

AD 

106 

40 

IF (QMAX) 33,33,41 

AD 

107 

41 

ARG4=SQRT ( SQRT ( RHAX/QMAX) )*ARG4 

AD 

108 


IF (HMAX-ARG4) 33,42,42 

AD 

109 

42 

SENS=0. 

AD 

110 


GO TO 34 

AD 

111 

43 

ARG4=HM IN 

AD 

112 


SENS=1 • 

AD 

113 


GO TO 3 

AD 

114 

44 

RETURN 

AD 

115 


END 

AD 

116- 
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O' U1 OJ 


SUBROUTINE DTLNEO C N« A , C,X , DET ) 
DIMENSION A(20,20>, C<20), X{20>, 
DET=0 • 0 
DO 2 1 = 1 1 N 

cl ( ii-cii * 

DO 1 J=1,N 

1 AL( I , J )=A ( I , J ) 

2 CONTINUE 
IR=G 
PRD^l.O 
ALT=I .0 
DO 15 J=I , N 
ALT A=- ALT 
BIG=0.0 
DO 8 1*1, N 
I RH 0= 1 

IF IIR) 21,6,4 
I RHO= IRHO+1 
IF (I-MCIRHOn 5,8,5 
IF (IR-IRHO! 21,6,3 
ALT A=-ALTA 
V = ABSC AL( I,J) I 
IF (V-BIG) 8,8,7 

7 BIG=V 
ALT=ALT A 
M(IR+U=I 

8 CONTINUE 

IF (BIG I 21,21,9 

9 I=M ( IR+1 1 
PRD-PRD*AL C I , J ) 

CL(I)-CL(n/AL(If J> 

IF (J-N) 10,16,21 

10 L*J + 1 

DO 11 K=L , N 

11 AL(I,K)=ALU,K)/AL(I,J) 

I R= IR+1 

DO 14 I U= 1 * N 

DO 12 IRHO-1 , IR 

IF (IU-M(IRHGI) 12,14,1 2 

12 CONTINUE 

DO 13 K=L, N 

13 AL( I U, KIGALI IU,K)-ALUU,J)*ALU,K) 


AE 1 

AL ( 20,20 ) , CL ( 20) , MC20I AE 2 

AE 3 
AE 4 
AE 5 
AE 6 
AE 7 
AE 8 
AE 9 
AE 10 
AE 11 
AE 12 
AE 13 
AE 14 
AE 15 
AE 16 
AE 17 
AE 18 
AE 19 
AE 20 
AE 21 
AE 22 
AE 23 
AE 24 
AE 25 
AE 26 
AE 27 
AE 28 
AE 29 
AE 30 
AE 31 
AE 32 
AE 33 
AE 34 
AE 35 
AE 36 
AE 37 
AE 38 
AE 39 
AE 40 
AE 41 
A.E 42 
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CL< IU)»CLUU)-ALC IU,J»*CL( I ) 

AE 

43 

14 

CONTINUE 

AE 

44 

15 

CONTINUE 

AE 

45 


GO TO 21 

AE 

46 

16 

DET=ALT*PRD 

AE 

47 


I=M(N) 

AE 

48 


X(N)=CL( I I 

AE 

49 

IT 

IF (J-1I 20*21tl8 

AE 

50 

18 

L=J 

AE 

51 


J=J-1 

AE 

52 


I=M ( J) 

AE 

53 


xu)=cnn 

AE 

54 


00 19 K=L,N 

AE 

55 

19 

X( JI=X( JI-AL< I,KI*X(K) 

AE 

56 


GO TO 17 

AE 

57 

20 

DET=0. 0 

AE 

58 

21 

RETURN 

AE 

59 


END 

AE 

60- 
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o o 


SUBROUTINE HTP (ITSVO 


AF 

1 



AF 

2 



AF 

3 

COMMON /AAA/ A( 6 , 6 ) , A12 (16 , 12 ) , A22 ( 7, 7 ) , ALP ( 6) , ALPT( 1 6 , 

5) , BMT ( 16 ) , 

AF 

4 

1C 1 6,16 ) ,C HI 16,2), CP 116) ,DEGI!6) ,H(16) , JAT( 16 , 5) , JPH ( 1 6 ) 

,K AT( 6 ) ,K0D 

AF 

5 

2E(16),RA(16,2 ) , RBOt 16, 2 ) , RC ( 16, 2 ) , RD( 16,2 ) , RD1( 16,2 ) , RE < 16, 2 ) ,RE1 < 

AF 

6 

316, 2), SOI 16) ,TB(3) ,TU( 16,2) »TU2 (16,2), VN( 17) , VNE ( 16 ) , VNT( 17 ) , VNU ( 1 

AF 

7 

46,6>,W26(6) ,W3(16),Y{16),RF(16,2),RC1 1 16,2),HS, IG, IGMS, 

IGMSP 9 IGP f I 

AF 

8 

50N, IS, ISP, ISPNGP,ISPNG2,N,NG,NP, PRESS, W27, RHO, WM, SYU, ISPNG 

AF 

9 

COMMON /RAD/ YY 1 100 ) , TEE ( 100 ) , FHV ( 20 ) , NHV, NY,C2, 1 Y 


AF 

10 

COMMON /RAO/ XNN(7,100) 


AF 

11 

COMMON /RAC/ N IHVC, FHVC ( 50 ) , AHV ( 50 ) , AHVLt 20 ) 


AF 

12 

COMMON /RAD/ Cl,C3,C4,FLG,C5,FLGl 


AF 

13 

COMMON /RAO/ YDELT, DELT0,FL1,FL2 


AF 

14 

COMMON /RAD/ GEEI 8) , EP S (8 ) , NU( 20 ) , ND( 70 ) ,FF ( 70) , GAMP { 70 ) , WOL ( 20) , F 

AF 

15 

1HVM (20 ) , FHVPI 20 ) 


AF 

16 

COMMON /CIONCL/ F ( 1 00, 1 0) , F2 ( 100, 1 0 ) , HVL( 70 ) ♦ EP,K2 , K1 , IFL, IYCON, IQ 

AF 

17 

1L,WMI,BIJ( 100,10) , GMIN ( 100, 10) ,GPLU( 100,10) ,IAEO 


AF 

18 

DIMENSION NICN(IOO) 


AF 

19 

DIMENSION PRES(IOO) 


AF 

20 

DIMENSION HH(IOO) 


AF 

21 

ALP ( 1 )*6.9375 


AF 

22 

ALP (2 ) =5. 5 


AF 

23 

ALP (3 )= 1 .4375 


AF 

24 

I TSW =2 


AF 

25 

.1 ON=l 


AF 

26 

I S=3 


AF 

27 

N =6 


AF 

28 

M26 ( 1 ) =0. 54862E-03 


AF 

29 

W26(2)=14. 007451 


AF 

30 

W26(3) =15.999451 


AF 

31 

KAT(i)=0 


AF 

32 

KAT ( 2 ) =7 


AF 

33 

KAT(3)=8 


AF 

34 

DEG I ( 1 )=— 1 .0 


AF 

35 

DEGI( 2 )=1 .0 


AF 

36 

DEG I (3 ) = 1 ,0 


AF 

37 

DO 1 J= 1 , 6 


AF 

38 

KODE( J)=2 


AF 

39 

DO 2 J= 1 , 6 


AF 

40 

DO 2 K= 1 , 5 


AF 

41 

ALPT( J,K)=0.0 


AF 

42 
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2 


JAT ( J* K)=0 

AF 

43 

ALPTtl, 11=1.0 

AF 

44 

ALPT( 2 » 1 )=2 .0 

AF 

45 

ALPTl 2 » 2) =2 .0 

AF 

46 

JAT (2»2)=7 

AF 

47 

ALPT ( 3 i 1 ) =1 .0 

AF 

48 

J AT (3*11 = 8 

AF 

49 

ALPT( 4* 1 ) =1.0 

AF 

50 

ALPT ( 4 , 2 ) = 1 .0 

AF 

51 

JAT (4*2 )=8 

AF 

52 

ALPT (5*11=1 *0 

AF 

53 

J AT( 5* 1 1=7 

AF 

54 

ALPT { 6* 1 ) =1 .0 

AF 

55 

ALPT(6,2>=1.0 

AF 

56 

JAT(6,2)=7 

AF 

57 

RA( 1* l ) =0 • 0 

AF 

58 

RC( 1*1 )=3. 9483 

AF 

59 

RD( 1 * 1 ) =6 .46872E-04 

AF 

60 

RE( 1, l)=-1.57124E+06 

AF 

61 

TU( 1,1) =5000. 

AF 

62 

RAl 1*2 )=0 .0 

AF 

63 

RC(1,2>=7. 71015 

AF 

64 

RD(1*2)=1. 42221 E- 04 

AF 

65 

RE( 1* 2 ) =— 3 .59273E+07 

AF 

66 

TU( 1,2)=24000. 

AF 

67 

RCl ( 1*1 ) = 6 .4483 

AF 

68 

RD1(1,1)=6.46 87 IE— 04 

AF 

69 

RE1 ( 1* 1 )=-1.57124E+06 

AF 

70 

TU2 (1* 1 )=5000 • 

AF 

71 

RCK 1,2)= 10. 2102 

AF 

72 

RD1(1»2)=1.4222 IE— 04 

AF 

73 

RE1 ( 1* 2 )=-3 .59274E+07 

AF 

74 

TU2( 1 » 2 ) =24000. 

AF 

75 

RA( 2, 1 )=0 .0 

AF 

76 

RC(2,1)=24.9173 

AF 

77 

RD(2,1)=1 .08123E-03 

AF 

78 

RE(2,1)=-2.21 576E+06 

AF 

79 

TU( 2* 1 )=5000. 

AF 

80 

RA( 2,2 ) =0.0 

AF 

81 

RC(2,2>=30»6229 

AF 

82 

RD(2,2)=3 .01129E-04 

AF 

83 

RE( 2*2 )=— 5. 19808E+07 

AF 

84 

TU( 2*2 )=24000. 

AF 

85 
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RCi ( 2 » 1 >=28.2374 

A F 

86 

ROl ( 2 , 1 )=1.32 899E-Q3 

A F 

87 

RE1 ( 2 , 1 >=— 2. 1 492E+06 

AF 

88 

TU2 (2,1 >=4000. 

AF 

89 

RC1(2,2>=33.5174 

AF 

90 

RD1 (2, 2 )=4.39119E— 04 

AF 

91 

RE1(2»2) =-3. 1 4737 E+07 

AF 

92 

TU2 (2, 2 >=24000. 

AF 

93 

RA(3»1)=3. 73033E+05 

AF 

94 

RC(3,1)=20.0611 

AF 

95 

RD( 3»1 ) =6 .473Q7E— 04 

AF 

96 

RE ( 3, 1 ) =— 1 . 57026E+06 

AF 

97 

TU( 3,1 >=5000. 

AF 

98 

RA(3»2>=3. 73033E+05 

AF 

99 

RC(3,2>=23.5272 

AF 

100 

R0(3,2)=l .7478E— 04 

AF 

101 

RE(3,2)=-3.23627E+07 

AF 

102 

TU( 3,2 >=24000. 

AF 

103 

RCI (3, 1>=22.7653 

AF 

104 

ROl (3,1 >=5.88991E-04 

AF 

105 

RE1 (3, 1 >=- 1.72 548 E+06 

AF 

106 

TU2(3,1 >=6000. 

AF 

107 

RC1(3,2)=26.3126 

AF 

108 

RD1(3,2>=1 .93706E— 04 

AF 

109 

RE 1(3, 2) = — 5. 03986 E+07 

AF 

110 

TU2 ( 3, 2 > = 24000 . 

AF 

111 

RA(4,l)=0.58985E+05 

AF 

112 

RC(4,1 >=20.8169 

AF 

113 

R0( 4, 1 )=0.655447E-03 

AF 

114 

RE ( 4, 1 ) =— 0. 16 2744 E+07 

AF 

115 

TU( 4, 1 ) =5000. 

AF 

116 

RA( 4, 2 ) =0 • 58985E+05 

AF 

117 

RC(4,2)=24.5103 

AF 

118 

RD( 4,2 )=0 . 158277E— 03 

AF 

119 

RE( 4, 2 ) =—0.351471 E+08 

AF 

120 

TU( 4,2 > =24000 . 

AF 

121 

RCI (4,1 >=23.3381 

AF 

122 

ROl (4, 1 )=0.659562E— 03 

AF 

123 

RE1 (4,1 >=-0.1 5558 7E+07 

AF 

124 

TU2 (4, 1 >=5000 . 

AF 

125 

RCI (4,2 > = 26.8167 

AF 

126 

ROl (4,2 > = 0. 186999E-03 

AF 

127 

RE1 (4, 2 )=— 0.322302E+08 

AF 

128 
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TU2(4,2)=24000. 

AF 

129 

RA( 5,1 )=4.48051E+05 

AF 

130 

RC(5,1)=20.609 

AF 

131 

RDC 5,1) =6. 56764E-04 

AF 

132 

RE ( 5, 1 ) =— 1 . 6422 6E +06 

AF 

133 

TU(5,1)=5000. 

AF 

134 

RA( 5,2 )=4 .48051E+05 

AF 

135 

RCC 5,2 ) =2 4 .346 

AF 

136 

ROC 5,2 ) = 1 .56145E-04 

AF 

137 

RE ( 5, 2 ) =— 3 . 5928E+07 

AF 

138 

TUC 5,2) =24000 • 

AF 

139 

RC1C5, 11=23, 3481 

AF 

140 

R01 (5, 1 )= 5, 99657E— 04 

AF 

141 

RE 1 ( 5 , 1 )=— 1, 70972E+06 

AF 

142 

TU2 C 5, 1 ) = 6000 . 

AF 

143 

RC1(5,2)=27.2144 

. AF 

144 

ROIC 5,2) = 1.53048E— 04 

AF 

145 

RE1 ( 5, 2 )=— 4.98219E+07 

AF 

146 

TU2 (5,2 ) =24000, 

AF 

147 

RAC 6, 1 ) =0 , 1125E+06 

AF 

148 

RCC 6,1)=19,8515 

AF 

149 

RDC 6, 1 )=0 .650822E-03 

AF 

150 

RE(6,1 )=— 0.156336E+07 

AF 

151 

TUC 6,1 )=5000. 

AF 

152 

R AC 6, 2 )=0 . 1 125E+06 

AF 

153 

RC(6,2)=23.2779 

AF 

154 

ROC 6,2 )=0.190986E— 03 

AF 

155 

REC 6,2 )=-0. 331003 E+08 

AF 

156 

TUC 6 ,2 )=24000 , 

AF 

157 

RC1 C6, 1 )= 2 2, 463 3 

AF 

158 

RD1C6, l)=0.622485E-03 

AF 

159 

RE1C6, 1 )=— 0.16534 5E+0 7 

AF 

160 

TU2 (6,1 )=6000 , 

AF 

161 

RC1C6, 21=26.0356 

AF 

162 

RD1 (6,2 )=0. 229541 E-03 

AF 

163 

RE1 (6, 2 )=— 0.5 0169 6E+G8 

AF 

164 

TU2 (6, 2 )= 24000. 

AF 

165 

RETURN 

AF 

166 

END 

AF 

167- 
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SUBROUTINE TRAP1 IN , X, Y ,C, ANS ) 

AG 

I 

DIMENSION X(l) , YU) 

AG 

2 

ANS=0. 

AG 

3 

DO 1 1 = 2, N 

AG 

4 

ox=xm-x(i-i) 

AG 

5 

ANS=ANS+.5*DX*(Y( IS)+Y( 1-1) ) 

AG 

6 

ANS=C*ANS 

AG 

7 

RETURN 

AG 

8 

END 

AG 

9 
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FUNCTION FLUT (A,B*C) 

AH 

1 

XN=A/B 

AH 

2 

R=IC/*00129)*29. 

AH 

3 

FLUT=2 ® 69 E19*XN*R 

AH 

4 

RETURN 

AH 

5 

END 

AH 

6 
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qo-4 N) h- on 


SUBROUTINE INOCAL (DX,X,S,T,L) A[ 

DIMENSION DXCI) * X(l) AI 

AI 

AI 

COMMON /AAA/ AC 6, 6 ) , A12 < 16 , 12 ) , A22 < 7 , 7 ) , ALPC 6 ) , ALPT C 16 , 5) , BMT < 16 ) , AI 
1C(6,16),CH< 16,2),CP(16) ,0EGI(6),H( 16) ,JAT{ 16,5) ,JP4<16) ,KAT<6) ,KQD AI 
2E(16) f RA{16 f 2),RB0(16 f 2),RC(16,2),RD( 16 , 2 ) , RD1 ( 16 , 2 ) , RE ( 1 6 , 2 ) , RE 1 ( AI 
316,2) ,SD( 16) ,TB(3) , TUI 16,2) ,TU2( 16,2) , VN< 17) ,VNE( 16) , VNT( 17) ,VNUC 1 A I 
46 ,6 ) ,W26{ 6 ) , W3 (16)*Y(16)*RF<16,2) , RC1 ( 16,2) , HS * IG, IGMS, IGMSP,IGP,I AI 


50N , IS , I SP , I SPNGP, ISPNG2 ,N , NG, NP, PRESS , W27 , RHO , WM» SYU, I SPNG A I 

COMMON /RAD/ YY ( 100 ) , TE E ( 1 00 I , FHV< 20 ) , NHV, NY, C 2 , I Y AI 

COMMON /RAD/ XNN(7,100) AI 

COMMON /RAD/ N IHVC, FHVC < 50 ) , AHV t 50 ) , AH VL< 20 ) AI 

COMMON /RAD/ Cl ,C 3, C4, FLG, C5 , FLG1 AI 

COMMON /RAO/ YDELT, DELTD, FL 1 , FL2 AI 

COMMON /RAO/ GEE( 8 ) , EPS ( 8 ) , NU( 20 ) , NDC 70 I , FF ( 70 ) , GAMP < 70 ) , WOL ( 20) , F AI 
1HVMI20) ,FHVP( 20) AI 

COMMON /CIONCL/ F< 100, 10) , F2 < 100, 1 0 1 , HVL( 70 ), EP , K2, K1 , I FL, I YCON, IQ AI 
11 *WMI f BIJ( 100,10) ,GMIN 1 100 ,10) , GPLUI 100,10 ) , I AED AI 

DIMENSION NICNUOO) AI 

DIMENSION PRES ( 100 ) AI 

DIMENSION HHUOO) AI 

IF (S) 10,1,9 AI 

SUM=0 • AI 

DO 5 I=K1,K2 AI 

J=I-K1+1 AI 

IF (F ( I Q1 , J ) * EQ *0 * ) GO TO 5 AI 

GO TO (2,3) , I AED AI 

E P“ GM IN(IQ1,J) AI 

GO TO 4 AI 

EP-GPLU ( I Q1 , J ) AI 

CONTINUE AI 

E = 1 « / ( ( T-HVLL II )**2+EP*EP) AI 

E^E*1.17E-29 AI 

SUM=SUM+E*(F2 ( IQ1 , J >-F2 ( I Y, J))*FLOAT< IFD/3.14159 A I 

CONTINUE AI 

IF (SUM-88.0) 6,6,7 AI 

DX= 1 *0— EXP ( SUM ) AI 

GO TO 8 AI 

DX=1.0 AI 

W MI = X AI 

RETURN AI 

9 CONTINUE AI 

10 RETURN AI 

END AI 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 
46- 
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SUBROUTINE TRAP ( N« X* Y f C, ANS ) 
DIMENSION ANSdl, XC1I, Yd) 

ANS { I )=0 • 

DO 1 1=2, N 

ox=xc n-xu-i i 

ANS C 1 1 =ANS( I-iI+.5*DX*m I)+Y(I-1))*C 

RETURN 

END 




SUBROUTINE MAS8AL i ¥N? I B,C, ALP, IRHQ, BD? IS,N, PR, IG t 

AK 

1 


DIMENSION VNC 17 ) , IBf 13), CC6,16), ALP(6), AC20,20) , BC20), XI20) 

AK 

2 


NP= N+l 

AK 

3 


ISP=IS+1 

AK 

4 


DO 2 LS = 1 , IS 

AK 

5 


BILSl-VNCNPI^ALPl LS ) 

AK 

6 


J = 1 

AK 

7 


DO 2 I = 1 , N 

AK 

8 


Bit S > = 8 ( LS !-C CLSt I)*VN II) 

AK 

9 


IF U-IBIJI) 2,1,2 

AK 

10 

1 

A(LS,JI=C{LS,I! 

AK 

11 


J=J*1 

AK 

12 

2 

CONTINUE 

AK 

13 


DO 3 J=1,ISP 

AK 

14 

3 

A C I SP t J )-0#0 

AK 

15 


DO 4 LS=1 » I S 

AK 

16 

4 

AILS, ISP)=-ALPI LS 1 

AK 

17 


BU SP ) = PR 

AK 

18 


DO 5 1=1, IG 

AK 

19 

5 

BUSPJ=8USP!-VN(I> 

AK 

20 


J=0 

AK 

21 


DO 7 1*1, IS 

AK 

22 


IF CIB(I)-IG) 6,6,7 

AK 

23 

6 

J=J*l 

AK 

24 

7 

CONTINUE 

AK 

25 


DO 8 1*1, J 

AK 

26 

8 

AUSPiINUO 

AK 

27 


CALL DTLNEQ ( ISP, A, B,X, V) 

AK 

28 


J=1 

AK 

29 


G*1.0 

AK 

30 


DO 13 1*1, NP 

AK 

31 


IF II-IBU)) 13,9,13 

AK 

32 

9 

V=BD*VN« I ) 

AK 

33 


W*ABS(X(JI) 

AK 

34 


IF CW-VI 12,12,10 

AK 

35 

10 

IF fV/W-G) 11,12,12 

AK 

36 

11 

G=V/W 

AK 

37 

12 

J=J+1 

AK 

38 

13 

CONTINUE 

AK 

39 


J=1 

AK 

40 


DO 15 1*1, NP 

AK 

41 


IF U-IB(J)) 15,14,15 

AK 

42 

14 

VNI II*VNC IH6*XiJJ 

AK 

43 


J=J+1 

AK 

44 

15 

CONTINUE 

AK 

45 


RETURN 

AK 

46 


END 

AK 

47- 


108 



USAGE 


Program Information and Deck Configuration 

The program D1250 was written in the FORTRAN IV language for the Control 
Data 6600 computer system under the Scope 3.2 operating system. The program requires 
a field length of 110 000 octal locations; the computing time depends primarily on the 
accuracy of the initial 6 0 and the number of iterations desired, KETEST. The following 
sketch shows the deck configuration needed for execution: 







For the sample case (U^, = 15.24 x 10^ cm/sec (50 000 ft/sec), and altitude of 60.96 km 
(200 000 ft), and KETEST = 12), the run time was 18.21 minutes. The results of a simi- 
lar case are discussed in reference 3. 

Input Description 

The input is loaded by using the FORTRAN IV NAMELIST. The input symbols are 
as follows with typical values given in parentheses: 

$EQRAD 

ALT altitude for AT62 
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AO VERB 
CCI 

cn 

CIMAX 

DELTX 

DUMBTIM 

EELE1 

EELE2 

EELT 

EITR1 

EITR2 

EMUREF 

EPSR 

EPS21 
I ALT 

IGEO 


ratio of semiminor axis to semimajor axis of an ellipse 
initial computing interval for INTI (0.05) 

initial linear integration step away from singular point of stagnation 
streamline; used in FIRST (0.05) 

absolute value of maximum computing interval in INTI 

an increment used in ITR1 so that is a reasonable 

DELTX 

approximation to derivative of f(x) (-1.E-3) 

initial guess for standoff distance = (0.78 for sphere) 

RHOl 

array of relative error for dependent variables used by INTI (l.E-3) 

array of relative zero for dependent variables used by INTI (l.E-4) 

array of NT values to return to the program from INTI at specific values 
(NT is computed in CONT) 

relative error criterion for ITR1 (l.E-4) 
absolute error criterion for ITR1 (l.E-4) 
molecular weight of cold air, g/mole (28.967) 
accuracy criterion for Rj iteration in MAIN (0.01) 

accuracy criterion used in PRO PIT for RHO iteration (l.E-3) 

option for atmosphere routine: 

= 1 call AT62 and input ALT 
* 1 input PIN, RHOIN 

test in GEO: 

1 
2 
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3 


sphere geometry 
ellipsoid geometry 
hyperboloid geometry 



ITEXT 


ITTEST 

KETEST 

KKKK 

PEEP 

PIN 
QI1 
QI2 
Q01 
Q02 ^ 

R 

RB 

RHOIN 

RHON 

R1 

SSPEC 

TEAN 


time history print option for INTI; if equal to 0, no printout is requested; 
if equal to 1, a printout is requested 

iteration limit in PROPIT (500) 

number of iterations desired on DELTN (10) 

radiation option: 

0 do not compute radiation 

1 adiabatic radiation 

2 nonadiabatic radiation 

coefficient in equation describing pressure variation in stagnation region 
(see ref. 3, eq. (36)) (-2.5 for sphere) 

free -stream pressure, dyn/cm^ 


curve -fit coefficients to radiation distribution around body, used for all 
but the final 6 0 iteration (see discussion on p. 42) 


gas constant, ergs/mole -°K 

body radius at axis of symmetry x = 0, cm 

free-stream density, g/cm^ 

initial guess for RHO, nondimensionalized 

initial guess for weighted heat flux, nondimensionalized (if value not 
known, use R1 = 0) 

print control option for INTI 

initial guess for T', °K 
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TCG 


accuracy criterion for iteration in RANH (0.001) 


UIN free -stream velocity, cm/sec 

XMAX limiting value in x-direction 

$ 


The input for the radiation portion of the program was supplied with the RATRAP 
code. Part of the data is permanent and, therefore, is not normally modified. The input 
is loaded by using FORTRAN IV NAMELIST. The input symbols are as follows: 


$RAD1 

GEE 

EPS 

FHVM 

FHVP 

FHV 

WOL 

NU 

ND 

HVL 

FF 

GAMP 


permanent data supplied with code for occupation number, density 
calculation, and line absorption coefficient 


NIHVC total number of continuum spectral points 
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FHVC 

AHV 

AHVL 

FL1 

FL2 

FLG 

FLG1 

NIC 

NICN 

IOPT 


NY 


TB 


spectral points for continuum calculation 
transmission factors for continuum calculation 
transmission factors for line calculation 

1 intensity calculation 

2 flux calculation 

1 intensity calculation 

2 flux calculation 

0 normal output 

± 0 diagnostic output 

0 intensity calculation 

* 0 flux calculation 

number of y/5 points where transport evaluated 

an array of NIC values of index location of y/5 points where transport 
is evaluated 

option in subprogram FEMP: 

1 temperatures supplied 

2 enthalpies supplied 

number of values of YY, PRES, TEE, and HH (maximum is 100) 
initial temperature estimate 


A listing of the NAMELIST input data appears in the output for the sample case. 


Discussion of Output 

The output of program D1250 consists of printing only. The NAMELIST input 
EQRAD is printed from MAIN and the NAMELIST RADI from PROD. If PIN and RHOIN 
are found by the atmosphere routine, they are printed in MAIN. Also printed in MAIN 
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are the initial DELTA for each 6 0 iteration and the maximum x from the previous 
iterations. The remainder of the output is printed only during the final 5 0 iteration. 

The results from the integration routine are printed by CONT. The frequency of the out- 
put in CONT is controlled by a print option in INTI. If SSPEC is zero, control is returned 
to the program after every acceptable integration step. The distribution across the shock 
layer is written by PROD and the radiation results are written by RADFLUX. The output 
symbols are as follows: 

AT62: 

RHOIN p 

r oo 

PIN p 

OO 

DELTN 6 0 

X x 

Stagnation -point conditions: 

UO 
RHOO 
AO 

RHOOAO 
HO 
PO 

DWDX 

Shock standoff distance: 

DELTD dimensional 5, cm 

DELTA/RB s/r b 


u 0 

P 0 

a 0 

p o a o 

Hn 


^0 

dw 

dx 
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Distribution across shock layer: 


ETA 

V 



PDUM 

P 



STAENT 

h 



CAPH 

H 



AUDUM 

u 

velocity gradient when 

PROD is called by MAIN 


u 

velocity when PROD is 

called by CONT 

V 

V 




Radiation output: 

PATH LENGTH 

PRESSURE 

TEMPERATURE 

ENTHALPY 

DENSITY 

DENSITY (ND) 

QRYPC body and shock continuum radiation heat flux, watt s/cm ^ 

QRYPL body and shock line radiation heat flux, watts/cm^ 

QRYO total radiation heat flux at body, dimensionless 

QRY1 total radiation heat flux at shock, dimensionless 

Properties at first step off stagnation line: 

X x 
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DELTA 


6 


W w 

10(2) to 10(5) Ij )0 

SMALRB r b 

Body and shock properties: 


(j 


X 

DELTA 

W 

10(2) to 10(5) 
SMALRB 
11(2) to 11(5) 
RBX 
THETAB 


x 

6 

w 

r b 

Rg at any x-coordinate 

% 


0 


(i 


Q 

DDELDX 

DWDX 

DI0DX(2) to DIODX(5) 

DSMALRBDX 

PRATIO 

QRATIO 


Q 

d6 

dx 

dco 

dx 

dx 

5b 

dx 

PO 

PSTAG 


ratio of heat flux at body to stagnation heat flux 


(j : 


2, 3,4, 5) 


2, 3,4, 5) 


2, 3,4, 5) 


2, 3,4, 5) 
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Output for Sample Case 
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QU * —0*302 E+0 1 * 



T0+3Si86T'0- = TOO 
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EPS = 0*0, 0.2384E+01, 0. 3576E+01, 0.10452E*02, 0.11877E+02 

0*1 3002E+02 1 0.0* 0.0, 


*> 
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BODY AND SHOCK PROPERTIES 
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RHO IS SMALLER THAN ALLOWED 



Error Diagnostics 

The following table gives the error diagnostics as the message and the corre- 
sponding subprogram which prints the message: 


Message 

Diagnosis 

Subprogram 

ICODE =, iteration did not converge 

Nonfatal 

MAIN 

Iteration limit in RANH program stopped 

Fatal 

RANH 

RHO is smaller than allowed 

Nonfatal 

f THEP 
V THER 

RHO is larger than allowed 

r Fatal 
Nonfatal 

THEP 

THER 

H is zero or negative 

Nonfatal 

f THEP 
THER 

H is larger than allowed 

Nonfatal 

( THEP 
V THER 

Iteration limit reached in PROPIT 

Nonfatai 

PROPIT 

FEMP BLEW 

Fatal 

RADFLUX 

N-R did not converge TBJA = 

Nonfatal 

FEMP 

T did not converge TGUS =, DIG =, LK = 

Nonfatal 

FEMP 


Langley Research Center, 

National Aeronautics and Space Administration, 
Hampton, Va., July 13, 1971. 
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SUBROUTINE INTI 


Language : FORTRAN 

Purpose : To solve a set of ordinary differential equations. 

Use : CALL INTI (II, N, NT, Cl, SPE C, CIMAX, IERR, VAR, CUVAR,DE R, ELE 1, ELE 2, EL T, 
ERRVAL,DERSUB, CHSUB,ITEXT) 

n INTI is a closed subroutine composed of an initialization section and 

an integration section. The user is required to enter the initialization 
section before he starts his first integration step. The above calling 
sequence is used for both initialization and integration, with the value 
of the code word n determining which of the two sections of INTI will 
be entered. 

The user must set n = 0 in order to initialize. 

During initialization the derivatives will be evaluated by using the ini- 
tial values of the variables, but no integration will occur and control 
will be returned to the calling program. When INTI is called with 
n > 0, entry is made to the integration section. Upon each entry to 
INTI, the subroutine stores a 1 in II so that the user need not supply 
a value of II > 0 for repetitive integration. 

Besides serving as a means for specifying the entry point to INTI from 
the calling program, II can also be set to specified values in CHSUB: 

2 The user will store the integer 2 in II if the answers in CHSUB 
are not acceptable to him and he wishes to recompute the 
answers using a shorter interval. This shorter interval must 
be stored by the user in Cl. It must be smaller than the com- 
puting interval just used. 

3 The user will store the integer 3 in n if he wishes to return to 
the calling program. The answers for the interval are con- 
sidered acceptable to the user and will be transferred to the 
VAR array (explained below) by INTI. 

In DERSUB, II may be set to: 

4 The user will store the integer 4 in II if he wishes to discontinue 
calculation of the present interval and return to the calling 
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N 

NT 

Cl 


SPEC 


CIMAX 
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program. On return to the calling program, the answers at the 
beginning of the interval will still be in the VAR array. 

If the user does not set II to a value in either CHSUB or DERSUB, 

II will always be 1 upon the return to the calling program. 

An integer value supplied by the user which is the number of differen- 
tial equations to be solved. INTI is compiled to solve a maximum of 
20 equations but may be recompiled for larger values of N if necessary. 

An integer value supplied by the user which is the number of values in 
the ELT block described later in this writeup. INTI is compiled with 
a maximum of 10 values in the ELT block but may be recompiled for 
more values if necessary. 

A floating-point value supplied by the user which is the computing 
interval that INTI will use initially. Cl must be a signed value - posi- 
tive if integrating forward, negative if integrating backward. Upon 
entry to CHSUB, Cl will contain the computing interval that INTI will 
use for the next step unless it has to take a short interval to hit an ELT 
value or a SPEC value described below. The computing interval used 
on the present step is available in CHSUB as the algebraic difference 
between CUVAR (1) and VAR (1). Since the subroutine is used on a 
binary computer and the interval variation is a halving and doubling 
process, Cl should be a power of 2. 

A floating-point value supplied by the user which specifies how often he 
wishes INTI to return control to the calling program so that the user 
may print his results. 

SPEC = 0.0: Control will be returned after every acceptable integration 
step. 

SPEC > 0.0: SPEC is the absolute value of the specified increment of 
the independent variable for which the user desires control returned 
to the calling program. 

The first printout is made at the initial value of the independent vari- 
able. The next return is at the nonzero integer multiple of SPEC 
closest to the initial value of the independent variable. The remaining 
returns occur at values which have been updated from this point by the 
increment given in SPEC. The return times generated by the incre- 
ment given in SPEC are not altered by an intervening return due to an 
ELT value (explained subsequently in "ELT"). 

A floating-point value supplied by the user which is the absolute value 
of the maximum computing interval that will be used. This value will 
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CUVAR 


be used if the doubling process would extend the computing interval to 
a value larger than CIMAX. CIMAX should be set to 0.0 if there is no 
desired maximum. 

An integer value supplied by INTI as an error code. It must be checked 
at every return to the calling program. It may have the following 
values: 

1 A normal return, no error. 

2 The ELT block is not monotonic in the direction of integration. 

3 The variables have failed to meet the local-truncation-error 
requirements nine consecutive times. The answers at the begin- 
ning of the interval are still in the VAR array. 

4 The variables have failed to meet the local -truncation-error 
requirements at least nine times over the last three intervals. 

An acceptable answer has been reached, however, and is in the 
VAR array. v 

A double -precision one -dimensional array containing the independent 
variable followed by the N dependent variables. The user must store 
the N+l initial values (in the double -precision mode) in the array for 
initialization. INTI will store the new values of the variables in VAR 
after each integration step when they are accepted by the user in CHSUB. 
The elements of the VAR block can be printed out in the calling pro- 
gram in accordance with the user's specification in SPEC. 

A double -precision one -dimensional array which is given values by 
INTI for two purposes. INTI will store in the same order as the VAR 
array the values of the independent variable and N dependent variables 
at which the derivatives are to be evaluated in the DERSUB subroutine. 
Although CUVAR must be in a double -precision array in INTI to main- 
tain the "partial double -precision mode" of computation, the evaluation 
of the derivatives should be in single precision. Two suggested ways 
of doing this are as follows: (1) Consider CUVAR as a single -precision 
array of 2(N+1) elements in the DERSUB subroutine and, when using 
the ith element in a computation, assign to it the subscript value (2i-l), 
and (2) at the beginning of the DERSUB subroutine, transfer from 
CUVAR to some newly defined single -precision array and evaluate the 
derivatives using the latter. 

INTI will also store the tentative answers after each integration in the 
CUVAR array before calling CHSUB so that the user can check these 
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DER 


ELE1 


ELE2 


ELT 


ERRVAL 


values to decide to accept or reject the answers. If accepted, the 
CUVAR values will then be transferred to the VAR array. The decision 
as to whether the computation in the CHSUB subroutine should be done 
in single or double precision is a function of the individual application. 

In most cases single precision is adequate and can be accomplished by 
applying suggestions (1) and (2) in the preceding paragraph to the VAR 
and CUVAR arrays. 

No values need to be initially stored in CUVAR. 

An N+l single -precision array in which the user will store the deriva- 
tives evaluated in DERSUB. The derivatives should be arranged by the 
user in DERSUB in the same order as the VAR block so that DER (2) 
will be the derivative of the variable stored in VAR (2), for example. 
DER (1) will be unused. The derivatives must be computed using val- 
ues of the variables which have been stored in CUVAR (not VAR) by 
INTI. To avoid unnecessary double -precision computation, the user 
should apply the suggestion described under CUVAR. 

A one -dimensional array of N values supplied by the user each of which 
is the upper bound of local relative truncation error for the respective 
dependent variables. If the error for any variable exceeds its respec- 
tive ELE1 value, the computing interval is halved and the integration 
restarted at the beginning of the present interval. If the error for all 
the variables is less than 1/128 of their respective ELE1 values, the 
computing interval is doubled for the next integration step. 

A one -dimensional array of N values supplied by the user which repre- 
sents a small value or "relative zero" for the respective dependent 
variables. If the absolute value of any of the variables is less than its 
respective ELE2 value, the relative error criterion for that variable 
will not be applied. 

A one -dimensional array of NT values supplied by the user which are 
values of the independent variable at which the user specifically desires 
control returned to his program. The values in the ELT block must be 
monotonic in the direction of integration, or an error return will be 
given by INTI. 

A one -dimensional array of N elements in which INTI stores an esti- 
mate of the local truncation error for each of the N dependent variables. 
The relative errors are computed from these values and compared with 
the specified ELE1 values. 
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DERSUB 


CHSUB 


ITEXT 


A subroutine written by the user which will be called by INTI to eval- 
uate the derivatives. The derivatives must be stored in the DER array. 
INTI will call DERSUB to evaluate the derivatives with the values of 
the variable it has stored in the CUVAR array. 

These evaluations should be done with the use of single -precision arith- 
metic. The name given to the DERSUB subroutine must appear in an 
EXTERNAL statement in the calling program. The user may return 
to the calling program by storing a 4 in II. 

A subroutine written by the user to allow certain logical control. After 
each integration step, INTI will make available to the user in CHSUB 
the tentative answers in the CUVAR array. The VAR array will con- 
tain the last accepted answer (that is, the value of the variables at the 
beginning of the interval). Whenever the user specifies that the answers 
are acceptable, the values in the CUVAR block are transferred to the 
VAR block. In CHSUB the DER block will contain the values of the 
derivatives evaluated with the present CUVAR block. The user has 
three options: 

(1) Not change II. In INTI, II = 1 denotes that the user has accepted 
the answers in the CUVAR block. Upon entry to CHSUB from 
INTI, II always equals 1. 

(2) Set II = 2. The user does not accept the answers and wishes to 
recompute the interval using a new computing interval which he 
stores in Cl. This computing interval must be smaller than the 
computing interval just used. This new value of Cl will now be 
stored by INTI as the normal computing interval for the sub- 
sequent integration steps. 

(3) Set H = 3. The user accepts the answer but wishes to denote a 
condition that he can test in the calling program. Control will 
be returned to the calling program with the answers in the 
CUVAR array transferred to the VAR array. 

The name given to the CHSUB subroutine must appear in an 

EXTERNAL statement in the calling program. 

An integer code word supplied by the user which gives him the option 
to have INTI print out a time history of the computing interval and the 
reasons for its variation. This printout should be requested only for 
problems which must be rerun because of unsatisfactory results the 
first time. 
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ITEXT = 0: No printout is requested. 

ITEXT = 1: A printout is requested. 

Restrictions : See arguments listed under CALL statement. 

Method : INTI is a fifth-order integration subroutine. The classical fourth-order 

Runge-Kutta formula is applied in conjunction with Richardson's "extrapolation to the 
limit" theory. INTI is a variable -interval -size routine in which the interval is varied 
to meet a specified local relative truncation error. 

Accuracy : The variable -interval -size mode of logic is used to make available an esti- 
mate of the local relative truncation error which is then controlled as explained in 
the ELE1 block discussion. 

Roundoff error is controlled by use of the "partial double -precision mode of computa- 
tion" as explained in reference (a) of this subroutine. 

Reference : (a) Henrici, Peter: Discrete Variable Methods in Ordinary Differential 
Equations. 

Storage : 2703g locations 

Subroutine date : August 1, 


John Wiley & Sons, Inc., c.1962. 


1968 
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SUBROUTINE ITR1 


Language : FORTRAN 

Purpose : To solve the single equation of the form x = f (x) for one real root by the 
Newton -Raphson iteration method. 

Use : CALL ITR1 (X,DELTX, FOFX, E 1,E2, MAXI, ICODE) 

An initial guess supplied by the user. On a normal return to the calling 
program from ITR1, X contains the root. 

An increment supplied by the user so that + & - is a 

reasonable approximation to the derivative of f(x). 

A function subprogram to evaluate f (x). 

Relative error criterion. 

Absolute error criterion. 

A maximum iteration count supplied by the user. 

An integer supplied by ITR1 as an error code. This code should be 
tested by the user on return to the calling program. 

ICODE = 0: Normal return. 

ICODE = 1: Maximum iteration exceeded. 

ICODE = 2: Derivative = 0. 

Restrictions: A function subprogram with a single argument x must be written by the 
user to evaluate f(x). The name given to the FOFX subprogram must appear in an 
EXTERNAL statement in the calling program. 

Method : The Newton-Raphson iteration technique (ref. (a) of this subroutine) is used 
where 

X n+ l = qX n + (1 - q)f(X n ) 


V W-l(» n-l) 

X n - X n _i 


X 

DELTX 

FOFX 

El 

E2 

MAXI 

ICODE 
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Accuracy: The iteration process is continued until either of two convergence criteria 
are satisfied. These criteria are 

(1) if f(X n ) ^ ep then 

(2) if |f(Xn)| < e p then |f(X n ) - X n s e 2 

Reference : (a) Scarborough, James B.: Numerical Mathematical Analysis. Fourth ed., 
Johns Hopkins Press, 1958, p. 192. 

Storage : 137g locations 

Subroutine date: August 1, 1968 
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SUBROUTINE AT62 


Language : FORTRAN 

Purpose : AT62 approximates the U.S. Standard Atmosphere, 1962 (ref. (a) of this 
subroutine). Computes density in slugs /ft pressure in lb/ft^, temperature in 
degrees Kelvin, and velocity of sound in ft/sec at any geometric altitude z in the 
range -16 500 feet % z ^ 2 320 000 feet. 

Use : CALL AT62 (Z,ANS) 

Z Geometric altitude in feet. 

ANS A one -dimensional array that contains the results. 

ANS(l) Density in slugs/ft^ 

ANS(2) Pressure in lb/ft^ 

ANS(3) Temperature in degrees Kelvin 
ANS(4) Velocity of sound in ft/sec 

Restrictions : Range: For altitudes below -16 500 feet the values of density, pressure, 
temperature, and velocity of sound are not valid. The concept of the velocity of 
sound in the atmosphere becomes essentially meaningless at altitudes in excess of 
300 000 feet. To point out this limitation, the velocity of sound at altitudes above 
300 000 feet is set equal to the velocity of sound at 300 000 feet. For altitudes above 
2 320 000 feet, density, pressure, and temperature are set equal to their respective 
values at 2 320 000 feet. 

Method : The equations and techniques are identical to those used in computing the U.S. 
Standard Atmosphere, 1962 (ref. (a) of this subroutine). 

Accuracy : The tables in the referenced publication were computed with the IBM 7094 
electronic data processing system using some double -precision arithmetic. In con- 
verting the routine for the Control Data series 6000 computer systems all double - 
precision arithmetic was eliminated. Accordingly, there may be slight differences 
between the results of the converted subroutine and the tables. 

Reference : (a) Anon.: U.S. Standard Atmosphere, 1962. NASA, U.S. Air Force, and 
U.S. Weather Bur., Dec. 1962. 

Storage : 1654g locations 


Subroutine date : August 1, 1968 
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